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PREFACE 


This document reports a study carried out by Aerophysics 
Research Corporation during the period from April to August, 

1969, under Contract NAS 2-5383. The National Aeronautics and 
Space Administration Technical Monitor was Mr. Richard H. 
Petersen, Mission Analysis Division, Moffett Field, California. 
Mr. D. S. Hague functioned as Aerophysics Research Corporation's 
Project Leader for the study. 

The trajectory program and basic steepest-descent optimi- 
zation program which applies to multi-staged vehicles whose 
stage points occur at fixed times were originally developed 
under U.S. Air Force Funding Contracts AF 33 (616) -6848 and AF 
33 (657)-8829. Mr. B. R. Benson of the Air Force Flight Dynamics 
Laboratory sponsored these developments in the period from 
1959 to 1964. The optimization program was further extended 
under NASA Contract NAS 2-3691. Mr. Hubert Drake of the Mission 
Analysis Division, Moffett Field, California monitored the study. 

Major extensions to program capability were undertaken in 
the period from 1964 to 1965 by the author while at McDonnell- 
Douglas Corporation, St. Louis. These extensions included 
solution of problems involving simultaneous determination of 
both optimal stage points and time varying control histories, 
the multiple-arc problem, and the extension to two-vehicle 
problems when the maneuvers of one vehicle are pre-determined, 
the "maneuvering target" problem. The program delivered to 
Ames Research Center under the prese.“t study does not include 
the optimal staging and maneuvering ..arget capability. However, 
in the interest of information di ss - .nination, the present report 
collects in one source the analytic approach employed in the 
point mass trajectory equations, +h basic optimization formu- 
lation, the multiple-arc extension, and the maneuvering target 
extension. 

Past contributors to the progra a development include 
Mr. Robert L. Mobley, now with the I md Corporation, who pro- 
grammed the optimization formulation and 

Mr. Robert C. Browne, McDonnell-Dougias Corp. , St. Louis 

Mr. R. V. Brulle, McDonnell-Dougias Corp., St. Louis 

Mr. A. E. Combs, McDonnell-Dougias Corp., St. Louis 

Mr. K. N. Easley, now with Aerospace Corporation 

Mr. Ken Geib, McDonnell-Dougl as Corp., St. Louis 

Mr. G. D. Griffin, now with Aerospace Corporation 

Mr. G. G. Grose, McDonnell-Dougias Corp., St. Louis 

ur . H. L. Rozendaal, now with Lockheed Electronics, Houston 

Mr. F. W. Seubert, McDonnell-1 ’ougias Corp., St. Louis 

Mr. N. E. Usher, McDonnell-Dougias Corp , St. Louis 
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Mrs. Jane Yonke of Aerophysics Research Corporation prepared 
the present report. 

The steepest-descent atmospheric trajectory optimization 
programs described in this report are currently being extended 
to the general two-vehicle trajectory optimization problem. In 
the general case both vehicles exert time varying control with 
either cooperative or conflicting objectives. This work will 
be reported separately at a later date. Inquiries concerning 
the development should be directed to Mr. B. R. Benson of the 
Air Force Flight Dynamics Laboratory. 
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ATMOSPHERIC AND NEAR PLANET TRAJECTORY OPTIMIZATION 


BY THE VARIATIONAL STEEPEST-DESCENT METHOD 


by Donald S. Hague 
Aerophysics Research Corporation 


SUMMARY 


The variational steepest-descent method is described in 
detail. The analysis is first developed for a system consis- 
ting of a single-stage. Performance, constraint, and trajec- 
tory termination criteria may be any functions of terminal 
state and time. 

The analysis is then extended to multi-stage systems. 

Here, a stage is bounded by two points selected from initial 
conditions, points of state derivative discontinuity, and 
terminal conditions. In the multi-stage analysis the stage- 
times replace the independent variable time as the independent 
variable time itself thus becomes an additional state variable. 

The three-degree point-mass equations of motion for a 
vehicle maneuvering in the vicinity of a central planet are 
developed in rectangular rotating coordinates. Planetary 
characteristics include up to four harmonics in the gravita- 
tional field, flattening of the polar axis, layered atmosphere, 
and wind structure. Vehicle characteristics include generalized 
aerodynamic and thrust descriptions in terms of up to six 
control variables . 

Weighting matrix or control variable metric tensor defin- 
itions are presented in some detail. The problem of false 
convergence induced by ill-chosen weighting matrices is dis- 
cussed in terms of an order of magnitude analysis. 



INTRODUCTION 


Trajectory optimization by the Steepest-Descent Method is 
now a routine performance estimation at several government re- 
search establishments and major aerospace concerns. The computer 
program utilized for trajectory optimization studies in this 
report is capable of determining optimal three-dimensional 
flight paths for a wide variety of vehicles in the vicinity of 
a single planet. Atmospheric effects may be included, if 
desired. Past program applications include flight path opti- 
mization of 

a. High performance supersonic aircraft 

b. Spacecraft orbital transfer rendezvous and 
re-entry 

c. Multi-stage booster ascent trajectories 

d. Boost-glide re-entry vehicles 

e. Advanced hypersonic cruise aircraft 

f. Air-to-ground missiles. 

Optimal control can be. determined for any combination of 
the time varying variables 

a. Angle-of-attack (or pitch angle) 

b. Bank-angle 

c. Side-slip 

d. Throttle 

e. Two thrust orientation angles 

All the commonly employed terminal performance and constraint 
criteria may be specified. Inequality constraints may be im- 
posed along the vehicle flight path. 

Several options are available for specification of vehic 1 ~ 
aerodynamic and propulsive options. Data and vehicle charac 
teristics option can be modified at preselected stage points. 

An arbitrary number of stage points may be specified. 

Planetary characteristics are nominally set to those of 
the earth. Up to four gravitational harmonics may be specified. 
Nominal planetary atmosphere employed is the 1959 ARDC. A 
variety of wind specification options are available. An ellip- 
soidal planetary shape may be specified. 

The original trajectory optimization program is described 
in References 1 and 2. Equations of motion employed are 
described in References 3 and 4 . Some past applications are 
described in References 5 and 6. An extension of program capa- 
bility is described in Reference 7. An extension to '-imultan- 
eously determine both optimal time varying control and discrete 
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stage points together with some applications are described in 
References 8 and 9. A guidance and control application, the 
so-called lambda-guidance scheme, is reported in Reference 10. 

The optimization program of References 1 and 2 employs 
a second-order prediction scheme and several control variable 
"weighting matrix" options to assist convergence of the steepest- 
descent algorithm. These two features have also been included 
in a recently developed trajectory optimization. Reference 11. 
They are also retained as convergence options in an extended 
version of the program of Reference 1 and 2, reported in Ref- 
erence 12. 



THE STEEPEST DESCENT METHOD 


Problem Statement 

Point mass motion is governed by three second order 
differential equations of position together with a first 
order differential equation governing the mass. By suit- 
ably defining additional state variables, it is possible 
to reduce these equations to a set of first order differ- 
ential equations. Point mass motion is, therefore, governed 
by a set of first order differential equations. The form 
of these equations is 


= jf ^n(t), a B (t), tjj 


n - 1,2 N 

■ *= 1,2 M U) 


That is, there are N state variables whose derivatives 
x n (t) are defined by N first order differential equations 
involving the state variables, together with M control vari- 
ables, a m (t) • and t» independent variable itself. 

Constraints may be imposed on a set of functions of the 
state variables and time at the end of the trajectory. In 
this case, a set of constraint functions of the form 

kj = j*p( x n (T) ' T )j * ° 

p- 1,2 P (2) 

can be constructed which the final trajectory must satisfy. 

Any one of the constraints may be used as a cut-off function 
which, when satisfied, will terminate a particular trajectory. 
The cut-off function can, therefore, be written in the form 

Q-a^T), Tj = 0 (3) 

and determines the trajectory termination time T. In all, 
then, when the cut-off function it. included, there are 
(P + 1) end constraints. 


Finally, it may be that some other function of the state 
viriables and time at the end of the trajectory is to be 
optimized. Hence, a pc, off function 





\ 


(4) 


which is to be maximized ^ 


ed, can be constructed. 
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Now, suppose that a nominal trajectory is available. The 
requirements of this trajectory are modest; it must satisfy 
the cut-off condition. Equation (3) , but it need not optimize 
the pay-off function or satisfy the constraint equations. To 
generate this nominal trajectory by integrating Equations (1) , 
the vehicle characteristics , the initial state variable values , 
and a nominal control variable history must be known. Once 
this nominal trajectory is available, the steepest descent 
process can be applied. To do this, the trajectory showing the 
greatest improvement in the pay-off function, while at the 
same time eliminating a given amount of the end point errois 
as measured by Equations (2) for a given size of control vari- 
able perturbation, is obtained by application of the Variational 
Calculus . 


Equations (2) provide an end point error measure, for they 
will only be satisfied if the end points have been achieved. 
Therefore, any non-zero ip p represents an end point error which 
must be corrected. A convenient measure of the control variable 
perturbation can be defined by the scalar quantity. 


DP 2 



[«a(t)j[v(t)Jj«a(t) 


dt 


a) 


where W is any arbitrary symmetric matrix. In the case where 
all control variables have a similar ability to affect the 
trajectory, W is taken equal to the unit matrix, and DP 2 be- 
comes the integrated square of the control variable perturbations 
6a (t). It might be noted that if Equation 5 is to have meaning, 
it is essential that all control variables have the same dime. - 
sions. To meet this condition, the control variables can be 
expressed in non-dimensional form. 

The constraint on control variable perturbation size repre- 
sented by Equation (5) is an essential element of the steepest 
descent process; for the optimum perturbation will be found by 
local linearization of the non-linear trajectory equations about 
the nominal path. To insure validity of the linearized approx- 
imation, the analysis must be limited to small control variable 
perturbations by means of Equation (5) which provides an integ- 
ral measure of the local perturbation magnitudes. 


Single Stage Analysis 

The steepest descent process has been outlined above. To 
implement this method, an analysis of all perturbations about 
the nominal trajectory must be undertaken. In the present 
report, all perturbations will be linearized; only first 
order perturbations in the control and state variables will be 
considered. The objective of the linearized analysis is 



determination of the optimum control variable perturbation in 
the sense discussed in the previous section. 


Denoting variables on the nominal trajectory by a bar 
<* B (t) | nominal * jot,u(t)j 

and 

jx^t/J nominal * | (7) 

where there are M control variables and N state variables. 

Now consider a small perturbation to the control variable 
history , 6a (t) ; this in turn will cause a small perturbation in 
the state variable history, 6x(t). The new values of the vari- 
ables become 


and 

| * |«(t)| + jia(t) | 

(8) 


► = I x(t)| + j6x(t) j 

(9) 

The nominal state variable and perturbed state 
histories can also be written as 

variable 

jx(t)j * jx(t 0 )| ■ 

t- Jf^x (t), a (t), t^j dt 

(10) 

| x (t)| = |x(t Q )j + ^ 4 ^ X> a + , t^| dt 

(ID 


Subtracting Equation (10) from Equation (11) and using 
Taylor's expansion to first order. 



M - {at)} - £ j |L . • «*” j « ■ {•«*>} 


(12) 


where 

f - t ^x(t), a (t), (13) 

and where the repeated index indicates a summation over ail 
possible values. Differentiation leads to 


d ' 

jt: 


ix(t) 


d? 


6x“ + 2L i«* 


6 


(14a) 



or in matrix form 

St { ax(t) } = [ F ] W + [ G ]{ 5a } 


(14b) 


where 




and G. . 
ID 



(15) 


Here the ( : , j )*"* 1 element lies in the i *"* 1 row and column of 

the matrices; F is an N x N matrix and G is an N x M matrix. 


The effect of these perturbations on pay-off, cut-off, and 
constraint functions must now be determined. A general method 
for obtaining these effects, known as the 'adjoint method,' 
Reference 13, is to define a new net of variables by the equations 

[x(t)]* - [f( t )] ’ [x ( C ) j (16) 

By specifying various boundary conditions on the X, the 
changes in all fv ictions of interest can be found in turn. To 
show this pre -multiply Equation (14) by X' and Equation (16) 
by 6x', transpose the second of these equations and sum with 
the first giving 

<•*>! * [*]’H • + WI°]|*-j 

-MM H U7) 

which may be written as 

j^(Mx)J-[l]'[o]H < 18 > 

Integrating Equation (18) over the trajectory 


- {x'sxKo ■ 


dt 


(19) 


Now def-'.ne three distinct sets of X functions by applying 
the following boundary conditions at t = T: 


{x (T)| 




(20a) 

(20b) 

(20c) 
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Equation (16) may now be integrated in the reverse direction 
(i.e., from T to t ) to obtain the functions, {A<+,(t)}, {Xjj(t)}, 
and (X^(t)J. 

Substituting each of these functions into Equation (19) in 
turn and noting that 


Lv t >JM ■ l£JH- **« 

(21a) 

|V T) JW - }- »B t , T 

(21b) 


(21c) 

follows that 



(22a) 

®t=r *Cl x aj[ G ]H dt + [ X Q (t o>J{* <(t o ) } 



(22b) 

Ht=T = ^J[ X ^[ G ]{ a “} dt + [V*o)] { 6x(t o)} 



(22c) 


Now, Equations (22) give the changes in pay-off function, 
cut-off function and constraint functions at the terminal time 
of the nominal trajectory; however, on the perturbed trajectory, 
the cut-off will usually occur at some perturbed time, T + t. T. 

In this case, the total change in the above quantities becomes 


d <t> 


dft 



5t N[°]M dt + [ x / t o ) J{* c(t «> ) } + * (T) AT 

it l Xn J[ G ]W dt + hMM + ^ (T) at 


T * * 

JtM[ G ]H dt + [ x * (to) ] 


(23a) 

(23b) 

(23c) 
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Equations (23) supply the change in pay-off, cut-off, and 
constraint functions on the perturbed trajectory. 


The time perturbation in Equations (23a) and (23c) may be 
elinunated by noting that, by definition of the cut-off function. 
Equation (23b) must be zero. 


iT - - + lwj{ 




(24) 


Substituting Equation C24) into Equations (23a) and (23c) 

* 

^ i • r -i / x . . - (25a) 


where 


|>oj [ G ]{*« + [x* n (t 0 )J |ix(t 0 j 

f W [ G ]W dt + { ix(t o>} 

o 

M-M 

M-M -*fe$ aL 


(25b) 


(26a) 

(26b) 


Equations (25) reveal the significance of the A functions, 
originally defined by Equations (16) and (20). At time t Q , 

gives the sensitivity of <{> (T) to small perturbations in the 
state variables at t Q . Similarly, X^q (t) measures the sensi- 
tivity of <f> (T) to small perturbations in the state variables at 
any time t. The sensitivity of the constraints dip to small 
state variable perturbations at any time is likewise defined 
by each row of the function X^^(t). 

A measure of the sensitivity of a trajectory to control 
variable perturbations can be obtained from the quantities Aa'^G 
and X^'^G. Consider a pulse control variable perturbation at 
time t', that is, 6(t-t'), where 6 is the Dirac delta function. 
With this type of control variable perturbation, it can be seen 
from Equations (25) that the changes in pay-off and constraint 
functions will be Aaq (t ' ) ' G (t ' ) and Xi|>fi (t ' ) ' G (t ' ) , respectively, 
for fixed initial conditions. 

In order to apply the steepest-descent process, the perfor- 
mance function change, Equation (20a), must be maximized; subject 
to specified changes in the constraints, Equation (25b); and a 



given size perturbation to the control variables. Equation (5). 
This can be achieved by constructing an augmented function in 
the manner of Lagrange which is to be maximized instead of d$. 
For the present problem, the augmented function is 


u \ HHH* + |w*»>J {*<*»)} 

+ti f T [aoj [w] dt 
Jt o 


(27) 


where the v are P undetermined Lagrangian multipliers, and y 
is a single undetermined Lagrangian multiplier. The objective 
now is to find that variation of the control variable history 
which will maximize U. 

Consider a variation of 6a , that is a 6 (6a) . Then, 
it is always possible to write any 6a distribution in the form 

1 6a|« | A(t)| k, or j^4aj = |A(t)J k (28) 

where A(t) prescribes the perturbation shape; and k, its mag- 
nitude. Now that part of Equation (27) which depends on 6a, 
the perturbation in the control variable, can be written in 
the form 

U . k ^M[ G ]{ A(t 4 “ + [ G ]» A<t) l at 

+ kt(V>j[w]H}« t29) 

J t o 

So that 

H - f L x +nJ [“]{*<*>} d ‘ [ x *n] [“]{**>} dt 

^0 


* [w] {*(t)} dt 


( 30 ) 
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or 


6U 



(31) 


where it has been noted from Equation (28) that 

6 (6a) = A ( t ) 6k (32) 

Now, since Equation (31) holds for any h(t) ± it follows 
that it is a general relationship. Further, for U to be an ex- 
tremal, 6U must be zero. 

If U has been_maximized by means of a control variable 
perturbation 6a, 6U must be stationary for all small pertur- 
bations to the 6a, that is, for all 6 (6a) . The only way in 
which Equation (31) can be zero fcr all 6 (6a) is for the 
coefficient of 6 (6a) to be identically zero. That this last 
statement is true follows from considering the case where, over 
some finite time interval between t Q and T, the coefficient of 
6 (6a) is, say, positive. If this were the case, we could 
choose a 6 (6a) distribution that was also positive in this 
same interval and zero elsewhere between t Q and T. It would 
follow that U was also positive, and, hence, U could not be 
maximum. A similar argument holds when 6 (6a) is negative ovej. 
any interval in t Q to T. Hence, the coefficient of 6(6 a) must 
be identically zero in the whole interval t Q <. t <. T. This 
argument is essentially based on that presented by Goldstein, 
Reference 14. It follows that 

|M* 


ii . 


(33) 



Transposing, noting that W is symmetric, and solving for 6a, 

(34) 


where 


— b[ w ] 1 M(M*M ' }} 

Substituting Equation (34) into Equation C25b ) 

V} + [ I w]{‘ , }| 

z o 

M- j* M'[°][ w ]' 1 [ o ]{ x J at 

“o 

For subsequent use define the integral 


and 


(35a) 


(35b) 


(36a) 


(36b) 


(36c) 


The multipliers v can be expressed in terms of the multipliers 
y by Equation (35a) 




(37) 


Substituting Equation (34) ini ■ Equation (5) 


>0X11^ LmI&C second term in the right hand 

^ + 2 MivM -J MM) 


(38) 


5?l 

Transposing the second term in the right hand side bracket 

(39) 

7 

Substituting Equation (37) in Equation (39) 
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and noting that is symmetrical gives 

v 2 !* 2 - i **-Lv*J[ i w]' 1 { i 4 t ‘ ta2 L M J[ I ^ 1 { dS } 1401 


So that 


On = L-fylpj 

'\dp 2T jd)5 j jd£J 


(41) 


Substituting Equation (41) into Equation (37), the remaining 
Lagrangian multipliers are obtained in the form 


M'-M " 1 


- ii^>j 1 i*»T 

DP 2 - Ld/3j [i^- 1 | d/3} 



(42) 


The optimum control perturbation is found by substituting 
Equations (41) and (42) back into Equation (34) and is 

{ 6ot } = l [w ] -1 [ G ]|{Vfi}'[Vfi] 1 M} 

v - / pP 2 - Ldfl l idfl{ 

+ [w]-l [g] [x^J [l„J* {dflj (43) 


With this equation the steepest-descent control pertur- 
bation has been determined. Perturbing the control variables 
according to Equation (43) gives the optimum change in the 
trajectory as discussed in the section entitled, "Problem 
Statement," with the added effect of changes in the initial 
value of the state variables included through the term in d£. 
The appropriate sign to use on the first term of equation (43) 
can be determined by evaluating d<j>. Substituting the optimum 
control perturbation into Equation (25a) results in the 
equation shown on the following page. 
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«♦ - y( r ** - m [ vF J w) ("* - m ta i h) 

* lvl[ I **]' 1 { <w } *l.WVj { (44) 

As the quantity in the radical must be positive to assure 
the change in <p is real, it follows that the negative sj'.n 
must be taken when minimizing the payoff function and the 
positive sign when maximizing the payoff function. 


An Alternative Analysis Using the Independent Variable 

for Cut-Off 

In the analysis of the previous section, it is implied 
that any function of the form 

fl(Xn(T), T) - 0 (45) 

will suffice to terminate the trajectory. While this is true 
in an analytic sense, in practice any function passing through 
zero more than once in the cut-off region may be difficult to 
employ for cut-off purposes. 

fl 4 



Figure 1. — Double Valued Cut-Off Function 


Figure 1 presents a nominal cut-off function history which 
decreases monotonically . The perturbed cut-off function 
history, shown dotted, behaves in a different manner in that 
it passes through zero twice in the cut-off region. As the 
trajectory must be integrated numerically , there is a danger 
that cut-off will occur the first time Q passes through zero 
ir.'tead of the second, thereby introducing both errors in 
the linearized perturbations and preventing the build-up of 
the anticipated cut-off function history. 



One method for overcoming this difficulty is to terminate 
the trajectory at the inearly predicted cut-off value of the 
independent variable, ’i + AT. This revises the analysis of 
the previous section; for by terminating the trajectory in 
this manner, a sma 1 ! error will exist in the value of the cut- 
off function, say Aft. In this case, Equation (24) must be 
modified to account for the cut-off function error. Allowing 
for this effect, we can write in place of Equation (24) 

( aii -fL [ x nJ [ c ] M dt 

6(t) 


(46) 


Substituting into Equations (23a) and (23b) 



and 


(47a) 

(47b) 


The additional term plays no part in the equations which 
result from taking a variation in 6qt;nence, the remaining 
equations ir this section up to Equation (34) are 'till cor- 
rect. Equation (34) mast now be substituted into Equation (47b) 
instead of Equation (25b) with the result 


M- 

Defining d(3* by 



(48) 


(49) 


and substituting d8* for d(3 in the remainder of the analysis 
of the previous section, the entire analysis remains correct 
up to and including Equation (43) . Equation (43) still gives 
the optimum control variable perturbation. The change in <j> , 
however, becomes 


- 15 - 
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MULTI-STAGE ANALYSIS 


Outline of Multi-Stage Analysis 

Trajectories in which some of the state variables or 
state variable derivatives have a discontinuity for some 
value, t', of the independent variable, t, are frequently 
encountered. Such a point will be called a stage point. 

That portion of a trajectory preceding a stage point is in 
a different "stage" to that following the stage point. The 
first stage will be that portion of a trajectory lying be- 
tween t Q and the first stage point; the s*^ stage will be 
that portion of the trajectory lying between the (s-1)*-* 1 and 
sth stage points. The final stage will be that portion of 
the trajectory lying between the last stage point and the 
final cut-off time. 


It is convenient in the analysis of multi-stage trajec- 
tories to define a new independent variable to replace t. 
This new independent variable, the stage time t, is defined 
separately for each stage in the following manner. Let the 
s tn stage commence at time t' . and terminate at time t' 
Then, s 1 s * 


T 


S 


= t 


fc s-l 


(51) 


so that when 


k - K-i - T s - 0 (52 > 

The s^ stage is terminated by a cut-off function G , 
assumed to be of the form 

n 8 - 0 8 ( X < T «>» T a) - 0 (53) 

where T g is the stage time at cut-off. 


The analysis of the preceding section, pages 4 to 16, 
no longer holds for a staged trajectory, unless the stage 
points are determined by cut-off functions of the form 


G (T ) = 0 (54) 

s s 

That is, the stages are of fixed length in the independent 
variable. 

For suppose the nominal trajectory has an s th stage 
lying in the region 


fc i-l 


<. t <. 


fc s 


(55) 
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1.U 

Then on the perturbed trajectory unless the s stage and the 
(s-l) th stage are terminated in the manner of Equation (54), 
the sth stage will occupy the region 

i s-i ' “;-i 1 1 1 ^ * 4t ; < 56 > 

by virtue of the perturbations in the state variables. This, 
in turn, will mean that estimates of the optimum a-perturbation 
based on the analysis of the preceding section, pages 4 to 16 
will be in error due to the fact that the F and G matrices 
are incorrect in the regions between 


(t s-l' 


t s-l 


At ;-i j 


and 


(t-, t; + At;) 

In such a situation, a new factor enters the optimization 
problem; for control may be exercised over the position of all 
or some of the stage points. In this case not just the opti- 
mum control variable perturbation is sought by the analyst but 
rather the optimum combination of control variable and stage 
point perturbations when considered simultaneously. The objec- 
tive in the following section is to obtain this optimum 
combined perturbation. 


Changes in Payoff and Constraint Functions 
in Combined Perturbation 

Given a nominal multi-stage trajectory, suppose the 
control variable histories and the stage point positions are 
simultaneously perturbed throughout the whole trajectory. 
Considering the first stage, the effect of this combined per- 
turbation, when the stage terminates, will appear as a modi- 
fication in the state variables values. This perturbation in 
the first stage final state variable values may be looked 
upon as a perturbation in the initial state variable values 
to the second stage. The combined effect of stage point, 
control variables, and initial state variable perturbations 
in the second stage will be to produce a perturbation in the 
state variable values at the second stage termination. These 
effects ca' ’e in like manner until the last stage is reached. 

Perturbatic in a typical stage are illustrated in Figure 2. 
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Figure 2. — Perturbations in the s*"* 1 Stage 


Consider the last stage. Since the result of trajectory 
perturbations ahead of this stage appear solely as initial 
state variable perturbations, the optimum combined perturbation 
along the whole trajectory can be found by optimizing the last 
stage in the presence of initial state variable perturbations 
which are a function of the previous stage perturbations. 

Consider the s*"* 1 stage of the perturbed trajectory. As 
noted above, the effect of perturbations in all stages pre- 
ceding the s^h appear as some perturbation in the initial state 
variable values of the s^* 1 stage, 6x . Suppose the s-* 1 stage 
is terminated on the nominal trajectory by some function of 
the form 

Then, provided the stage time, t s , is used as the independent 
variable instead of t, the change in any function of the state 
variables and t s at cut-off can be found by 

T i 

-f [ x * n s ] [ G ]{*«} dT 8 + [ x * n 8 (°)] { 5x 8 (0) } 


) 


(57) 


where ip Q is any function of the form 




( 58 ) 
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Kl-M-M ^ 

(59) 

Here Ai)> s and Aft s are obtained by integrating 
equations in stage time through the stage 

boundary conditions 

the adjoint 
subject to the 

[ X * S (T,)] " 

' a* si 

L* x jJ 

(60) 

n 

u 

c: 

/< 

\ asisj 
I axj j 

(61) 


So far, the cut-off function has not been perturbed; this can 
be achieved by terminating the trajectory when 


n + Afi = 0 
s s 


(62) 


instead of by Equation (53) . Perturbing the cut-off function 
will cause a change in the trajectory stage time at cut-off 
given by. 


A T s 



(63) 


The total change in i|» s at the termination of the s'" stage 
will then be given by 


/•Tj t | 

{**■(*8 + AT S >} «7 [ x *« b ] [ G ]{ 5a } dr s + [ x * fi 8 ( 0 )] / Sx 8 (°)} 


+ 



(64) 


Suppose the state variables are selected as \p s , so that 

{*.} - {*<T 8 )} (65) 

With this choice of i^ s . Equation (60) becomes 

KM ■ fe]’ •[*] (66) 

where I is the unit matrix. Denoting the Ai)i s resulting from 
this particular choice of boundary conditions by Ax s , it 
follows from Equation (64) that the change in state variables 
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at termination of the s th stage is 

| ix 8 (T g +AT S ) i ■/[■*•]'[ <>]{<«}dr 8 + [*^(0)] |«x 8 (0)} + {x(T 8 )JaT 8 
o 

(67) 

These perturbations are the state variable changes to the left 
of s + 1 stage point. 

The i x to the right of a stage point, 6x^ are not 
necessarily equal to those on the left, 6x'“', Figure 3, but 
a matrix P s can be defined which will transform the left-hand 
perturbations into the right hand ones. 

K:i} - [p.] W;i} ■ [*.] {»*,<*, 

( 68 ) 

Typically, for example, consider the case of a multi-stage 
boost vehicle with a fixed amount of fuel in each stage. The 
mass of the remaining portions of the vehicle at the commence- 
ment of the s^ stage is the sum of the empty weights of the 
remaining stages, together with the sum of the fuel contained 
in those stages. Perturbations in mass at the termination of 
the (s-l)th stage reflect changes in the burning time of that 
stage. It will usually be physically impossible to transfer 
any fuel remaining in the (s-Dth stage across the interface 
with the sth stage, and, hence, changes in the state variable 
of mass to the left of the stage point may fail to cause a 
corresponding change to the right of the stage point. In such 
a case, the P matrix will have a null row for that particular 
state variable. On the other hand, changes in the state vari- 
ables of position to the left of a stage point will always 
appear directly as changes to the right of a stage point. The 
co:. responding row in the P matrix will have unity on the diag- 
onal element and zero elsewhere. This is also true of the 
state variables of velocity, provided impulsive forces are 
absent at the stage point. 

Substituting Equation (67) into Equation (68) 

K4 - [ p 4f s M '[ 0 ]H a,s * M 0) ] 'K (0) } + 

-[*.]{{*.}♦ [W{ ,x - ( 0) MM 4T ‘) (69) 
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where 


r Ts 

N'J NWH"- 

O 


(70) 

= 

[w°>] 


(71) 

{**) = \ 



(72) 

|^* 4 +) 

< 



• 

X 

* 

* 8 

8 X 




s th stage* 


(s-l) th stage point 


s th stage point 


time 


th 


Figure 3. — Position of Functions Defined in the s Stage 


In some cases additional perturbations in the state 
variables to the right of a stage point may be specified. For 
example, returning to the case of a multi-stage booster, in a 
more sophisticated analysis of booster capability, one may 
wish to consider variations in the initial mass of fuel con- 
tained within each stage. Typically, in a given iteration the 
amount of fuel consumed in the s+lth stage may be either less 
than or greater than the total amount of fuel available in 
that stage. If this or a similar situation arises, the ini- 
tial amount of fuel contained in the stage must be adjusted 
in the next iteration. It is essential to have a mechanism 
within the optimization analysis which will permit these 
required changes to be specified. 

The P matrix, as described above, is unable to provide 
this mechanism, for the additional changes may clearly be 
functions of the state variable perturbations at the termin- 
ation of the (s+l)th stage rather than at the termination of 
the s^ stage. The P matrix is primarily introduced to con- 
vert changes in the state variables at the termination of the 
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s stage into state variable changes at the beginning of 
the fefl)*-* 1 stage. Accordingly, a set of additional state 
/ariable perturbations {Ax .} which are specified directly 
may be introduced. The complete expression for the state 
variable perturbations, therefore, becomes 

K4 ■ NjN + [x]’K <0) } + {‘ S K) + 


With Equation (73) the first objective of this analysis is 
achieved: a recursion formula which determines initial state 

variable perturbations in the (s+l)th stage when the per- 
turbations in the s^h stage are known and additional changes 
are directly specified in the initial values of the (s+l)th 
stage state variables. 


The recursion formula Equation (73) can be applied to 
each stage in turn, commencing with the first. In the case 
of the first stage, there are no perturbations in the initial 
state due to prior stages, but there may be perturbations to 
the state variables if a search for the optimum trajectory 
initial conditions is being made. 


These initial value perturbations will be some combi- 
nation of state variable vectors dictated by the particular 
problem under consideration. For example, suppose the op- 
timal launch point for a mobile mission is sought. The 
nominal launch point may be specified in terms of latitude 
and longitude. On successive iterations the launch point 
is perturbed towards the optimal position. The perturbation 
in latitude is a state variable vector having components 
corresponding to the position state variables and zeros 
elsewhere 


Here X , Y e 
the planet; 



/ d*_ 
<9X e 
d ± 
<9Y e 
d±_ 
dZ a 


\ 0 


(74a) 


and Z are rectangular coordinates rotating with 
as disSussed in a later section, pages 45 to 47. 


Similarly, 


the components for longitude changes are 



(?) 



(tI 


(74b) 
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In addition to this type of initial point perturbation, there 
may be changes resulting from previous iterations, Ax p . 
Combining both types of perturbation, the total change^in the 
first stage state variables is 

W} - { 4x i} ■ ( 75 ) 

Knowing the 6x^, the total initial value perturbations in 
the second stage may be computed from Equation ( 73 ) . 

{*x 2 } = [ p l]j{ K l} + [*xoJ {ax x } + {%}AT 1 | + {ax^ ( 76 ) 
Proceeding to the next stage, 

W * [ P 2]{{ K 2}*[ X *a,] {8*2} + {ix 3 } 

■ NH’NNMW 

* + [ P2 ] NT iTi 

+ N N]t Pp l NT M + H tm 'IN * {^^3} (77) 

and the next 

{ Jx 4 • [ P 3]|{ K 3} * IN] {**3} + { ^3} AT 3 } + { 4x 4 

= Q 3 Q2^lKi + + P3K3 

• • • 

+ x^AT^ + ^ 3 P 2 *2^2 + p 3 x 3^ r 3 

+ ^2*1 Ax l + ^3*2 Ax 2 + *3 Ax 3 + 

(76) 

where 

[*s]“ [ P# ][ Xxn s] (79) 

In general, the total state variable perturbation at the 
commencement of the (s+l)' :h stage is 

{sx 8+1 } = (Vs-1 *2 P 1 K 1 + *s*s-l ^ ♦ — 

+ Q s Pb-1^b- 1 + P s K s^ 

+ (ft S *8-l QgPiXjATi + Q b Q b . 1 <4 3 I > 2 x 2 AT 2 + 

+ Q S P 8 -1 X S -1 A ^S -1 + P S X S A ^ C ^ 

+ (^sVl— ft l Ax l + *s* 8 -l *2 Ax 2 + — - 

— + *8*8-1 Ajt 8-1 + *8 Ax 8 + AX 8+l) ( 80 ) 
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At the commencement of the last (Nth) stage 


w- & raw- &«".• tw 


( 81 ) 


where 


[*! 

1- 

^ N - 1%-2 • • • • 

^S+l^B > 

S < N-l 



= 


y 

S = N-l 

(82) 

{»! 

>■ 

%- l %-2 • • * * 

• 

, S < N-l 




• 

Vs 


, S = N-l 

(83) 

[<£ 

= 

%-l%-2 • ♦ • * 


, S < N 



=: 

I 


, S = N 

(84) 


Knowing the initial perturbations to the last stage. Equation 
(25) can be applied in stage time to find changes in the pay- 
off function <f>, and the constraints. 



o 



o 


(85) 


( 86 ) 


and on substituting for <$x N 

»♦ ■ %, *[ r H{s [ # 3M ♦ ’ j^^KS) <87) 

M - {*»#} ♦[*#>«] ’{i[ [*3 W * 5) KK * sMW) (88> 
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where 




l X *%J[ G ]{ ia } 


dr N 



(89) 


(90) 


It is convenient to jmbine the integrals, K s , through each 
stage into a single set of integrals throughout the complete 
trajectory. To accomplish this, define 


L^hJ[ A ®] [ Xx ^] ’ 


a < N 



Expanding Equation (87) 


(91) 

(92) 



( 94 ) 
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The first term, which is the summation of a set of integ- 
rals throughout each stage on the unperturbed trajectory, can 
be combined into one integral by reverting to the original 
independent variable t. That is 


**- / b-aJHM dt * 4T » * ^[ c "] W| 1951 

Similarly, 



(96) 


Equations (95) and (96) give the total change in paycff 
and constraint functions when the trajectory simultaneously 
undergoes perturbations in the control variable histories, 
stage point positions and initi al state variable values ir 
each stage. The sensitivity of payoff and constr; int functions 
to these variations is immediately apparent. For pulse 
variations in the control variables at time t = t’, the indi- 
vidual sensitivities of the payoff iunction are the elements 
of the row matrix 

j_ < ^ > 0t(‘ fc * )J = J 'J (97) 

The individual sensitivities of the constraint functions to 
control variable pulse perturbations are similarly the elements 
of the rectangular matrix 

[*«<**)] “ [*Vn] [ G ] os) 

The sensitivity of the payoff function with respect to 
stage point variations follows directly from the second term 
of Equation (95). If the s^h stage alone is perturbed by AT S , 




(99) 


Similarly, the constraint sensitivities with respect to 
star i point perturbations are obtained from Equation <9 6) 


{*AT S }“ [**■] W 


( 100 ) 
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Finally, the sensitivities to initial state variable 
value perturbations in each stage can be obtained from the 
last terms in Equations (95) and (96). The payoff function 
sensitivities for the s th stage are tie elements of 

|>.J- [°s] (10l > 

and the constraint function sensitivities for the s^h stage 
are the elements of the rectangular matrix 

[^xj " [°J] (102) 

In general two types of stage point must be considered: 
those whose perturbation is prescribed and those which are 
free to optimize. Hence, 



where the AT , are prescribed and the AT- are to be optimized. 
Substituting s Equation (103) into Equation (95) 

** ' /iMMM 4 ' 1 * MjsT 

* I'Hj.Sj® 2 *} * EKIMj uo4) 

Substituting Equation (103) into Equation (96) 



(105) 


where all the quantities specified directly are grouped to- 
gether in the term, 




(106) 


With Equations (104) and (105) expressions for the change 
in the payoff and constraint functions resulting from a gener- 
al perturbation of control variable histories and stage point 
perturbations are obtained. The optimal perturbation has yet 
to be found. This task is considered in the next section. 
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Derivation of Variational Equations 

The preceding section derived expressions for payoff 
and constraint function changes when a combined perturbation 
was introduced in the control variable histories, stage 
point positions, and initial state variable values in each 
stage. The perturbation which optimizes payoff function 
change while at the same time producing specified changes 
in the constraints {6ip} and the initial state variable 
values in each stage {AXg} will now be obtained. In 
order to obtain a meaningful solution, constraints are 
placed on perturbation magnitudes. Control variable per- 
turbations are limited in the manner of the preceding section 
by introducing the constraint 



(5) 


2 

A similar constraint (DT ) is introduced to limit total 
stage point perturbation. 


DT 


2 




) 


(107) 

Here the Vg are a set of weighting functions used to 
modulate optimal stage point perturbations. 

Proceeding as in the previous section , pages 10 to 13 
an augmented function is const' ucted in the manner of 
Lagrange and minimized (maximized) . In the present case, 
the augmented function is 
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T 



(108) 

where w is a Lagrangean Multiplier introduced for the stage 
point perturbation constraint. 

First, differentiate with respect to each stage point 
perturbation being optimized. This results in S equations. 


* L x# 0 *J{ b !} + 2w v » + W[v«n] {*£} 


(109) 


These expressions must disappear for U to be an extremal. 
Solving for the AT- 


AT-j 




( 110 ) 


Squaring both sides_of Equation (110) , multiplying throughout 
by V-, summing the S equations, and using Equation (107) 

S 
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2 2 ± 
E V- 4T/ ■ DT^ - — 


•dljN -wN^-H • wi<)|. 

(Ill) 

T . ranging Equation (111) and taking the summation into 
the matrix product 


• 2 - 


op( L ** * *l L wJ{*} * [«J [ L w]W) 


or 


where 


and 


wW - o } -|*|l w ]{ , } - 0 

L** " [ ®] { *"»} 

[ L +*J *| 3 "*J[ D ] 

M * 


( 112 ) 

(113) 

(114) 

(115) 

(116) 


Second, take a variation of 6a to Equation (108) in a similar 
manner to the preceding section 

»P JI ip 

*"<«) ■ t /[ A »nJ[°]{‘M dt * W/ * 2 | ‘i[H[“]{ 4 M dt * 


‘•[ A *«J[ G ] + [*J[ A *fi] [ G ] + 2 M [ a “J[ w ] * 0 


(117) 


(118) 


( 119 ) 
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Substituting into Equation (105) 

{ «•} --fejM-MW! 

where we define 


-J[ |>n|[ G ] [ w ] \ G ] W dt 

O 

t o 

[•w] * J t 


(120) 


( 121 ) 


( 122 ) 


(123) 


Transposing Equation (110) , substituting into Equation (120) , 
and using Equation (116) 

M- -5 s{M + W{»f-fe[ r ^]'[ D | r *«»]W 

Substituting Equations (114) and (115) into Equation (124) , 
multiplying throughout by 2wy and collecting terms on the 
left 

* 0 (125) 

Third, substitute Equation (119) into Equation (5) 

DrZ 't?l o [M N ][ °] M' 1 [°T{W * Mi* }| dt <i26> 

Using Equations (121) , (122) , and (123) , multiplying by 4y 2 
and collecting on the left 


* }• WMW ■ 0 (127) 
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Equations (112) , (125) , and (127) are the variational 
equations which must be solved to obtain the Lagrangean 
multipliers and, hence, the optimal perturbations. The 
solution to the optimal staging problem is immediately 
available in terms of the Lagrangian Multipliers. Equation 
(119) provides the optimal control variable perturbation 




(119) 


The optimum stage point perturbations are given by Equations 
(110). Rearranging and combining Equations (110) into one 
matrix equation 

* m h I v »J 1 f B 3 |{^%} + [v%] wj u28) 

where the columns of are the {B N } . The similarity 

of form between the expressions for {6a} and (AT S ) is 
immediately apparent. 

The optimal payoff function change is found by substi- 
tuting Equation (119) and (128) into Equation (104) 


1 ( 1 \ 

( 1 ( 1 1 ( .1 ) 

2/1 | + 



*Mj at„. * 

where the J terms are defined by Equations (121) to 
and the L terms by Equations (113) to (116). 

The equations presented here in this section constitute 
the general solution to the optimal staging problem. A 
closed form solution to these equations cannot be obtained. 

For straightforward elimination by Sylvester's Method leads 
to a high order polynomial equation indicating that a closed 
form solution is unattainable. Accordingly, in the next 
section a further assumption permitting attainment of a closed 
form solution is made. 


(129) 

(123) 


Solution Using Combined Step-Size Parameter 

The fundamental reason a closed form solution to the 
general optimal staging equations could not be obtained in 
the preceding section, pages 29 to 33 , was the introduction 
of separate Lagrangean Multipliers for the control variable 
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and stage point perturbations. If these perturbations are 
combined and a single step-size parameter defined by the 
expression 

f T r -i S 

DC 2 = J |te(t)J V(t) Ki«(t)}dt 4 £ (130) 

t o L L J Vxl 

this difficul .y can be avoided. It should be noted that the 
weighting functions, V-, must be dimensional quantities to 
insure compatibility between both portions of Equation (130) . 

Returning to the augmented function of Equation (108) 
it follows that in the last line to can now be replaced by yu.. 
This substitution can also be made in Equations (109) to (112) . 
Equations (112) and (127) can then be combined to obtain 

- (-W w - * [|>j * [ 1 **jj{* } 

(i3i> 

Substituting y for to into the equations leading to 
Equation (125) results in the expression 

fw} *{M} * M *[ L «] W ■ ° 

-* (132) 

It follows by comparison with Equations (39) and (35a) 
that y and {v} are given by Equations (40) and (42) provided 
the { dg } are replaced by { 6 r } and the I by the appropriate 
(J + L) . For example, 

► J** + 


The optimal control variable perturbation is then given 
by Equation (43) with the I replaced by the (J + L) and the 
X by the appropriate A . For example , 

► { A *n} 

Substituting y and {v} into Equation (128), the optimal 
stage point perturbations are obtained in the form 


34 




+ |" V ®J 1 [ B s] [voi] £w + 1 { 6r } ( 133 ) 

With these equations a general solution to the optimal 
staging problem has been obtained. The main difficulty in 
applying the results will undoubtedly lie in determining 
suitable weighting matrices. 
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POINT MASS TRAJECTORY EQUATIONS 


Basic State Variables 

Preceding portions of this report derived a successive 
approximation scheme for computation of optimum trajectories 
generated by a set of first order differential equations. The 
analysis is quite general and holds for trajectories generated 
by any set of first order differential equations. The object 
of the present section is to specialize the analysis to point 
mass vehicle trajectory problems. This will be accomplished 
when a suitable set of state variables, together with their 
derivatives, the control variables, and the forces associated 
with the control variables has been specified. First, a suit- 
able coordinate system is selected, and Newton's Laws in this 
system are used to define the vehicle's motion. 

Several suitable coordinate systems are available for 
point mass trajectory computations. The basic set of coordi- 
nates used in the present analysis is a rectangular set rota- 
ting with the earth, (X e , Y e , Z e ) . This coordinate system is 
illustrated in Figure 4. 



Figure 4. — Basic Coordinate System 
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The X e and Y e axes lie in the equatorial plane, the positive 
X e axis being initially chosen as the intersection of this plane 
with the vehicle longitudinal plane at t = t Q . Y e is 90° to 
the west of Xe, and Z e is positive through the South Pole. 
Denoting the radius vector from the center of the earth to the 
point mass vehicle by II , its magnitude is given by, 


l*|«Vx . 2 + Y . 2 + Z , 2 

(134) 


The angle between I and the North Pole is given by 

4> * - 90 - <fr L (135) 

where <p£ is the latitude of the vehicle. As a result of the 
earth's rotation, an observer in the (X e / Y e , Z e ) system would 
detect an apparent motion of the point mass if it were at rest 
in inertial space. In time At the apparent displacement of 
such a vehicle would be 

^apparent * R sin “p At (136) 

to the west. In vector notation 

*R«Iparent *fc*WpAt = -«p x RAt (137) 

This apparent displacement is independent of the vehicle's 
motion and exists whether or not the vehicle is at rest in 
inertial space. In general, then, to an observer in the rota- 
ting coordinate system, 


(ft R)e * ($R)ln*rtlal + ( $ R) apparent 
.*. (5 *) inertial “ (**)e +wp*^At 


(138) 

(139) 


Dividing Equation (139) by AT and taking the limit, it 
follows that 


((K^inertial " (<It)e 
or 

^inertial ' V e +<*> p xR 


(140 


(141) 


37 



The vector R in Equation (140) could equally well be taken 
as any vector; the arguments of Equations (134) to (141) still 
hold. Therefore, in general, for any vector quantity the oper- 
ational equality 

' (»). *“" X (142) 


can be defined. Applying Equation (142) to Equation (141) , the 
inertial acceleration is given by 


inertial 



(143) 


Now Newton's Law applies in inertial space so that in the 
rotating system 


l 

n 



+ *V Cw p x R 


(144) 


Here f is the total force acting on the vehicle. This vector 
equation can be expressed in component form using the relation- 
ships 


X e .| ♦ Y e .j + Vk 

(145) 

-Wp-k 

(146) 

V 1 + Fy e’J + f *e 

(147) 


Here l>j, end k are unit vectors aligned along the X e , Y e , and 
Z e axes, respectively. Equating components on either side of 
Equation (144) 
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Fv •• * 2 

- x e + 2wpY e - Wp* Xe 

(148) 

F v • * 2 

= Ye - a^Xe - Y, 

(149) 

P *e - z 


nr " * 

(150) 


These equations are not yet in a suitable form for the 
steepest descent analysis to be applied, for they are not in 
first order form. The transformation of Equations (148) to 
(150) into first order for:.> is immediately accomplished if 
the following quantifier are defined as state variables: 


*> Jt\ 




where 


(151) 


= u_.l + v«.i + v*.k 
cR- e * e 


(152) 


With this set of state variables the following expressions for 
the state variablt derivatives are obtained from Equations (148) 
to (152) 

*e “ «e 


L 


«e 

(153 

v e 

(154 

w e 

(155) 

ffe - 2«pv e + oJ Y* 

ID 

(156; 

F 

TT + 2w P u e + «p 2 Y e 

(157) 

F z 

z e 

» 

(158) 
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~e equations are in the same form as Equation (1) 
pic. the total forc r ‘ is a function of the state variables, 

a sj. i. control variaL . : j, and time. When the mass is vari- 
able, it too must be introduced as a atats variable. Any 
expression for *-he rate of change of mass of the form 

■ - i ^* n (t), a^(t), (159) 

may be used in the analysis. The above state variables, X e , 
Y e , u e , v e , w e , and m will be referred to as the basic 

state variables. In certain problems it becomes necessary to 
specify additional state variables; these will be discussed 
later in this section of the report. 


Control Variables 

The total force acting on the vehicle has three distinct 
sources; first, aerodynamic force as a result of interaction 
between the vehicle surfaces and the planetary atmosphere; 
second, gravitational force as a result of vehicle and planetary 
mass interaction; and finally, thrust forces introduced by the 
vehicle propulsion system. 

Befoie aerodynamic forces can be computed, the atmospheric 
properties, vehicle valocity relative to the atmosphere, and 
vehicle att must be specified. Atmospheric properties car. 

usually be specified as a function of altitude which, in turn, 
is a function of the state variables X e , Y e , Z^. vehicle velo- 
city relative to the atmosphere is also a function of the state 
variables, for u e , v e » and w e are the vehicle velocity compo- 
nents in a rotating system. The first and second factors 
determining aerodynamic forces are, therefore, f 'ctions of 
the basic state variables. 

The remaining factor entering into aerodynamic force 
determination, the vehicle attitude, is clearly not a function 
of the basic state variables. Fc - given the vehicle's position 
and velocity, we are still quite I. ee to specify its angular 
orientation in space. The angles which determine vehicle 
orientation may, therefore, be utilized as control variables by 
which aerodynamic forces may be modulated. Any set of three 
independent angles could be utilized for this purpose. Conven- 
tion suggests use of the vehicle angle-of-attack and angle-of- 
sideslip to orient the vehicle reference axis with respect to 
the velocity vector. Angle-of-attack, (a) , is the angle 

between the velocity vector and the vehicle reference axis when 
viewed in thj vehicle side elevation. That is in a rectangular 
coordinate system, x, y, z with x along the vehicle reference 
axis, pO"‘t’ve forward, y perpendicular to the vehicle plane of 
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symmetry, positive to starboard, and z completing a right hand 
system, a view normal to the x-z plane is considered. If u, v, 
w are the components of the vehicle velocity with respect to the 
atmosphere in this body axis system 

a - tan -1 (|j) (160) 



x 


V 


Figure 5. — Angle of Attack 


Sideslip angle (3) is the angle between the velocity 
vector and the reference axis when looking down on the vehicle 
planform, that is along the z axis. In this case, 

3 = tan -1 (£) (161) 



Figure 6. — Sideslip Angle 


Angle of attack and sideslip completely define th° attitude 
of tie vehicle with respect to the velocity vector. The third 
angle required to establish vehicle orientation in space is a 
rotation about the velocity vector. This last angle, bank angle 
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(B A ) , will be taken as zero when the vehicle plane of symmetry 
is vertical and the vehicle upright. Positive bank angle will 
be taken as a positive rotation about the velocity vector, as 
in Figure 7 . 



With the above set of angles to describe vehicle attitude, 
the velocity vector known and a given atmosphere, the aerodynamic 
forces can be completely specified. 

Returning to the second source of vehicle force, gravitation, 
from Newton's Laws, this is merely u. function of position and 
mass. It is, therefore, completely defined in terms of the state 
variables and, hence, introduces no new control variable. 

The final source of vehicle force, thrust from the propul- 
sion system, involves the atmospheric properties, either due to 
the atmospheric back pressure degrading the vacuum thrust or by 
virtue of the atmospheric fluid used in the combustion process 
which creates thrust. The propulsion unit efficiency may be 
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affected by Mach number and, hence, velocity so that thrust 
forces depend on the basic state variables of position and 
velocity in a similar manner to aerodynamic forces. If the 
propulsion system has a fixed orientation within the vehicle, 
the control variables introduced to describe aerodynamic forces 
suffice to describe thrust forces also. It may be, however, 
that the propulsion unit has a variable orientation within the 
vehicle. In this case, additional control variables to describe 
the relative position of the propulsion unit with respect to 
the vehicle are required. With vehicle attitude already spec- 
ified by a, 8 and B^, two additional angles are sufficient 
to orient the thrust. These may conveniently be taken as the 
cone angle from the reference axis, X-p, and the inclination 
about the reference axis, This latter angle will be 

measured positively aoout the reference axis and be zero when 
the thrust force is perpendicular to the port side of the veh- 
icle plane of symmetry, as illustrated in Figure 8. 



Figure 8. — Thrust Angles 

One other control variable for thrust remains to be speci- 
fied; this is the throttle setting, N, which serves to determine 
the propulsion unit power setting on variable thrust engines. 

In all then, to specify the forces acting on a point mass 
vehicle with a single propulsion unit, six control variables, 
a, 8, B^, At / <PTf and N,are required. If there is more than 

one independently controllable propulsion unit, additional X 
<j>T, and N must be defined. 
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Coordinates and Coordinate Transformations 


Local geocentric-horizon coordinates. — Components of 
the planet-referenced acceleration are integrated to obtain 
the planet-referenced velocity components (Xe, 2 e ). Vehicle 
position in this coordinate system is determined by integration 
of these velocities. Vehicle position ir. the planet-referenced 
spherical coordinate system will now be determined. The spher- 
ical coordinates are longitude, geocentric latitude, and dis- 
tance from the center of the planet. Angle "C" represents the 
change in vehicle longitude and may be written 

C “ 9 L C 0 L (162) 

Angle C is related to the vehicle position by the expression 

C = Tan -1 (~) (163) 

*e 

The relationships are illustrated in Figure 9. 



Figure 9. — Relation between Local-Geocentric, Inertial and 
Earth-Referenced Coordinates for Point-Mass Problems 


To describe body motion relative to the planet, a local- 
geocentric-horizon coordinate system is employed. The Zg-axis 
of this system is along a radial line passing through the body 
center of gravity and is positive toward the center planet. 

The Xg-axis of this system is normal to the Zg-axis and is posi- 
tive northward; Yg forms a right-handed system. Figure 9 shows 
the relation of this coordinate system to the other systems 
employed. 
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To locate the Xg-Yg-Zg axes with respect to the Xe - Y e -Ze axes, 
first rotate about Z e by an angle {180° + C) and then rotate 
about Yg through the angle (90° - <(,£) • The first rotation 
defines the intermediate coordinate system shown in Figure 10. 
The transformation is given by 


or 


V 


Cos (180° + C) Sin (180° + C) 0 


!x e 

l u 

S 

-Sin (180° + C) Cos (l80° + C) 0 





0 0 1 




ix’ 


—Cos C -Sin C 0 



% 

X 

Sin C -Cos C 0 


! Y e 



o 

o 




(164a) 


The second rotation is shown in Figure H. The transfor- 
mation matrix for the second rotation is given by 


1y 

A g 


Cos (90° - ♦l) 0 -Sin (90° - * L ) 


ix’ 


= 

0 10 


lv 

Y g 

l h 


Sin (90° - * L ) 0 Cos (90° - <p L ) 


1 7 
4 e 



-vi • , coordinate System TransforniaJion Figure n- Final Rotation in Transformation from 

- • ‘ ' rH •« Local-Geocentric Coordinates Earth-Referenced to Local-Geocentric Coordinates 
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In this analysis, a positive rotation is defined in the 
sense adopted for vector cross products in a right-handed system. 
That is, a positive rotation about the z-axis occurs when the 
x-axis rotates into the y-axis; positive rotation about the 
x-axis when y-axis rotates into the z-axis; and positive rota- 
tion about the y-axis when the z-axis rotates into the x-axis. 

The intermediate coordinate system (X 1 , Y , Z e ) is eliminated 
by the met od of successive rotation. The complete transfor- 
mation is given by 

Sin 0 -Cos -Cos C -Sin C 0 ly e 

0 10 Sin C -Cos CO ly 

e 

Cos <t> L 0 Sin <J> L 0 0 1 ±z e ( 165 ) 

I 

This can be reduced to the single transformation matrix 
-Sin Cos C -Sin Sin C -Cos <J>j\ lx e 

Sin C -Cos C 0 1 Y 

-Cos 4 >l Cos C -Cos <|>l Sin C Sin lz e (166) 

which defines a direction cosine set (l, j, k) by the equation 



lv 

A S 

iy 

6 




(167) 


Planet referenced velocity in the local-geocentric coordinate 
system is given by 


and 



(168) 

(169) 


Flight path angle* are computed by 
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and 


Y * sin 


-1 



(171) 


Here a is the heading angle and A is the flight path angle. 

Wind Axis Coordinates . — Aerodynamic and thrust forces 
for point-mass problems are c< -eniently summed in a wind-axis 
coordinate system, (X^, Ya, 2^ / . Since the equations of motion 
are solved in (X e , Y e , Z e ) coordinates, the wind-axis components 
of force must then be resolved into this basic system. 


When winds exist, defined by atmospheric velocity compo- 
nents along the local geocentric axes, vehicle velocity rela- 
tive to the atmosphere is the vector difference of vehicle 
geocentric velocity and wind velocity. The wind axis sysLem 
is then determined by the vehicle airspeed, V A , and the flight 
path angles relative to the atmosphere A a and °A* If wind 
velocity is zero, Va = Vg, Aa = A and a a = 0. If there is 
a wind, with velocity components (X gw , Y gw , Z gw ) , then 


V A «1 

1 ( v* 

f 8 j 

y a - 

. -1 

sin 

Q 

> 

II 

tan ^ 


* l VV 
[-<VV )/v a] 


e «v 


(y - i )/(x 

S gw g 


-X ) 
gw 


(172) 

(173) 

(174) 


Forces are first resolved from wind axes to the local 
geocentric coordinates. The wind axes are defined relative 
to the local geocentric axes by three angles: heading, 0 a? 

flight path attitude, X&, defined above; together with angle. 



Figure 12. — 

Relationship between Local-Geocentric Axes and Wind Axes 
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Appropriate transformations are 



X' cos o A sin o A 0 X g 

Y’ = -sino A cos a A 0 Y g 

-z- 0 0 1 z„ 

8 g 


( 175 ) 



(176) 


(177) 


The complete transformation from local geocentric horizon coor- 
dinates to wind axes then is 

X A cos y^cos <? A cos Y A sin ° A -sin y a 

Y a _ -sin a A c os B A cos a A cos B A cos Y A sin B A 

— + sin Y A cos c^sin B a + sin y. sin a^sin B A 

Z A sin o A sin B A -cos a A sin B A cos y a cos B a 

+ sin Ya cos 0 a cos B A + sin Y A s ^ n a A C0S B A 
which defines a direction cosine set 

x A r i s i t i x « 

Y a ■ r 2 s 2 t 2 Y g 

Z A r 3 S 3 t 3 Zg ( 179 ) 
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The resolution of forces from wind axes to local geocentric then 
becomes 


Fv 

A g 


r l r 2 r 3 


f x A 

% 

S 

S 1 ®2 »3 


F *A 

% 


t l t 3 


f Za 


(180) 


For the rotating-planet , the local geocentric components must 
be resolved into the X e -Y e -Z e system. The required direction 
cosines are given by Equation (168) 


P Xe 


i i h i 3 


Fy 

x g 

F *e 

S 

J]_ J2 J 3 


Fv 
1 g 

% 


k-L kg k^ 


F Z 


(181) 


The combined transformation from wind axes to local geocentric 
can be defined as a single matrix transformation 


P X e 


Oi 

X 

°2 

°3 


f x a 


”g x 

A e 

F *e 

= 

Pi 

P2 

*3 


f *a 

+ 

mg Y 

1 e 

% 
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92 

q 3 


f z a 


n «Z e 


(182) 


Body-axis coordinates. — Origin of this system is the ve- 
hicle center of gravity with x-axis along the geometric longi- 
tudinal axis of the body. Positive direction of the x-axis is 
from center of gravity to the front of the bod* . The y-axis is 
positive to starboard extending from the center of gravity in 
a water-line plane. The z-axis forms a right-handed orthogonal 
system. To permit the use of body (x, y, z) axes aerodynamic 
data, and to convert the body axes components of thrust to the 
wind axes system, a coordinate transformation must be made. The 
coordinate transformation shown in Figure 13 involves rotation 
first through angle of attack, a, then through an auxiliary 
angle, 6'. Noting that 


tan 8 ' * ^ cos a » tan 8 cos a 


(183) 
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the transformation is 

cos n 0 sin a x 

0 1 0 y 

-sin a 0 cos a z 

cos 6* sin 0' 0 

-sin 0' cos 0* 0 

0 0 1 

cos B f cos a sin 0* cos 0 ! r; 

-sin B f cor a cos 0 f -sin 0* 

-cin a 0 cos a 

which defines the (u, v, w) direction c 

X A U 2 Ug X 

Y a = v 1 v 2 v 3 y 

Z A w x w 2 w 3 z 

which define the force coefficient tram, formation 

-C D u i “2 J 3 ~ C A 

C Y/ = Vl u 2 u 3 C y 

” C L w l w 3 -C N (186) 

The relationship between body and wind-axes aerodynamic coef- 
ficients is now established. Note the negative directions of 
the coefficients relative to the axes. 

Inert ;. a.' c oordinates . — The selected inertial coordinates 
coincide w .r.h The earth references (X e , Y e , Z e ) system at time 
zero. At a la^er time they differ by the rotation of che earth. 
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Wpt. The transformation between inertial velocities and planet 
referenced velocities is derived as follows. 


Let R be the displacement of the point-mass, (See Figure 9). 
In inertial coordinates 


R = xl x + Y 1 y + Zl z 


(187) 


and 


(188) 


V = R = il x + Y1 y + zl z 
In planet-referenced coordinates 

R = X e I Xe + YeI Ye + Z e I Ze 

However, due to the rotation of the X e , Y e , Z e coordinate system, 
the velocity is 


where 


6R 

St 


H - + + 


(189) 


(1?0) 


The planet's rotation is about the Z-axis which it also the Z e - 
axis. Therefore, 


Wp - -ojpl z = -Wplz e 
and the required cross product is 


U) p X ^ = 


If Equations (188), (190), and (191) are substituted into 
Equation (189) , it follows that 

Xl x + Yly + zl z = (X e + WpY e )I Xe + (Y e - w p Xe)i Ye + <Z e )Iz e 

(192) 

The relation between the unit vectors in the inertial system 
and unit vectors in the planet referenced system are obtained 
by a single rotation about the Z-axis. 


x z c 


0 

x c 


0 

Y, 


- (Y e w p)ix e " ^X e Up)lY e 


(191) 
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( 193 ) 


The transformation from planet-referenced velocities to inertial 
velocities is made with the inverse of the matrix of Equation 
(193) and the component relations derived in Equation (192) 


• 

X 


Cos (Opt Sin (Opt 0 


Xe + “p^e 

Y 

= 

-Sin oipt Cos w^t 0 


Y e - (0pX e 

Z 


0 0 1 


z e 


(194) 


The components of inertial velocities are used to calculate the 
inertial speed of the body as 

V I * Vx 2 + Y 2 + Z 2 (195) 

Equation (195) is valid regardless of the inertial coordinate 
system involved. 


Local-geocentric to geodetic coordinates . — Positions on 
the planet are specified in terms of geodetic latitude and 
altitude (for a given longitude) while the motion of the body 
is computed in a planetocentric system which is independent of 
the surface. In the computer program, flight-path angle X 
and heading angle a are calculated with respect to the local 
geocentric coordinates. By definition Xq and op are angles 
measured with respect to the local geodetic. Although the 
maximum difference that can exist between the two coordinate 
systems is 11 minutes of arc, it may be desirable to know X^ 
and do more accurately than is obtained when measured from 
the local geocentric. 

It will be necessary to resolve the geocentric latitude 
to geodetic latitude for an accurate determination of position. 
Figure 14 presents the geometry required for describing the 
position of a point in a meridian plane of a planet shaped in 
the form of an oblate spheroid. 
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Figure 14. — Planet-Oblateness Effect on Latitude and Altitude 

It is apparent from this figure that the most significant 
difference between the geocentric referenced position and the 
geodetic position is the distance AB on the surface of the 
reference spheroid. The distance can be defined by a knowledge 
of the angle <j> L ; the geocentric latitude; <f>g, the geodetic 
latitude; the corresponding radii; and the distance OC. 

The relationship between the geocentric and geodetic lati- 
tude of a point on the surface of a planet which is an oblate 
spheroid is obtained as follows. The equation for the surface 
in a meridian plane is 



The tangent of the geodetic latitude can be found by determining 
the negative reciprocal of the slope of a tangent to this ellipse. 
The expression for this tangent is 

_ . 1 R e 2 Zb 

Tan *g = _ d£z) 1 ” rrr (197) 

dX B *p *B 

Note that Z R is a negative number in the northern hemisphere. 
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The tangent of the geocentric latitude of point B is 

Tan <p T = _ 

L s < 198 > 

Substituting Equation (198) into Equation (197) gives the 
required relation 

D 2 

Tan *g * __ Tan ♦l 

Rp 2 8 (199) 

The expression for the radius of the planet at point B in 
terms of the geocentric latitude of the point and the equato 
rial and polar radii is obtained by the rectangular to polar 
coordinate transformation 


-z B 


r *l Sin % 

g 6 




R 4>L CoS ^L a 
g 8 


( 200 ) 

( 201 ) 


and, solving for Ra by substituting Equations (200) and (201) 
into Equation (196)'"g gives 

— _ £s£p . — 

R * Lg = V R p 2 Si" 2 \ + Cos 2 ♦ “ 

( 202 ) 


cos <)> 


~ R e y/y[( R e/Rp)(tan <j) Lg /tan <J» L ) ] 2 sin 2 + cos 2 ^ 


COS tp 

It may be seen from Figure 14 that 


IPF' = OF' - OB' 


or 


h sin $ g = OP sin 4 >l - R^ g sin ^ 

Likewise 

FT" = OF" - 00" 


or 


(203) 

(204) 

(205) 


h cos 4g = OF cos * L - R* l cos 4 l (206) 

K * 

If Equation (204) is divided by Equation (206) and then the 
quotient is divided by tan <pr , there results 



tan <J> 

( — 

tan 


_ sin & 

[° p ( JITT > - ] / [«P < 


cor ’ 4 l , 

cos R<t>L ^ (207) 


or 


cos <j>L cos 

(R e /Rp) ( cos = (gin + [^ R e 2 - R p 2 )/ R p 2 ] [ R «J>L 4r /0P J ( 208 ) 


g 


Finally, if Equation (208) is multiplied by (Rp sin 4>L 9 )/ 
(Re sin (jjj,) , it follows that 9 

Rp « R^ sin 4 >t _ 

V If V = (^) *[l - (Hp/Re) ] ( r p - s ;„ ■ > W >(209) 


Let 


u » (R e tan <P L /Rp tan 4> L ) 
£> 


* (Rp tan <Pg/R e tan <f> L ) 

Then it follows from Equations (202) and (209) that 


u « & + [R /or] [U/ u 2 Sin 2 Z + COS 2 * L ] U - ( R / R e ) 2 l 
R 6 

e 

Equation (211) is solved by an iterative scheme. 

Then r u 

4> g = tan -1 [(-| — ) tan <*> L ] 


( 210 ) 


( 211 ) 


(212) 


The flight-path and heading angles corrected to the local 
geodetic latitude are computed by 




-(213) 


Since the magnitude of vector V g is equal to the magnitude of 

1 T * 


vector V, 
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and 


a D = Sin 


= Sin _P / -- 

\tTV + V/ 




V*g - +L>! * V/ 


( 214 ) 
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Auxiliary Computations 

In addition to the computations which can be made from the 
problem formulation as presented in preceding sections, several 
other quantities are available as optional calculations. 


a. Planet-surface referenced range, R D 

b. Great-circle range, R g 

c. Down- and cross-range, and 

d. Theoretical burnout velocity, V t j leo 

e. Velocity losses, Vp, V grav , V D , and Vml 

f. Orbital Variables and satellite 1 :.rget 

Planet-surfaced referenced range .?- The total distance 
traveled over the surface of the planet is computed as the 
integrated surface range. If the distance traveled by the 
vehicle over a given portion of the trajectory is 


% 



v g dt 


(215) 


then the curvilinear planet surface referenced r.mge is 


*■/ 5 * 

{ R 

t l 


Vg Cos y dt 


(216) 


The flight-path angle, X, is referenced to local geocentric 
coordinates for this computation. 


Great-circle range . — Great-circle distance from the launch 
point to the instantaneous vehicle position, R g , may also be 
required. Expressions for this distance are derived as follows. 


By spherical trigonometry, (see Figure 15 ) 


Cos » 


or simplifying 


Cos (90-4 > l )Cos( 90-4> Lo ) + Sin(90-4' L )Sin(90-^ Lo ) Cos(e L -e Lo ) ' 

(217) 


Cos » Sin Sin ♦ Lq + Cos ^ Cos ♦ Lq Cos (0 l - 8 Lo ) (2 18) 


Therefore, 

Rg ■ R' Cos ^ £sin ♦l Sin + Cos ^ Cos Cos (0 l“®L o ^J (219) 
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Figure 15. — Great-Circle Range 


However, since the planets are generally oblate spheroids, R' 
is not a constant radius. An approximation may be obtained by 
averaging the planet 1 s radius at the launch point and at the 
vehicle's position. Therefore, define the average radius, R' , 
as 


R’ 


R( i>L * R *Lq 
2 


( 220 ) 


and the surface-referenced great-circle range from the launch 
point to the vehicle is 


R g * 


R 4>L + R( ( , Lo 


Cos -1 r Sin<J>L Sin ^ + Cos ^ Cos ^ Cos(6 L -6 Lo )j 

L ° f 221) J 


Down- and cross-range . — Down- and cross-range from the 
initial great circle can be determined. The initial great 
circle is determined from the input quantities a Q , <f>L 0 , and 
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0L (see Figure 16) Then the cross range of a particular 
trijectory point is defined as the pe pendicular distance from 
the point to the initial great circle. The downrange is then 
the distance along the initial great circle from the initial 
point to the point P at which the cross range is measured. From 
the spherical triangle, Figure 16, the great circle range LF to 
the point F is computed by Equation (221) . 


The right spherical triangle LPF is solved for the down 
range, X D , and the cross range, Yq. 


Xn = B' 


COS 


■(; 


Cos LF 


Cos (sin -1 (sin LF sin ?)) 


) 


( 222 ) 


Y d = R' 


sin ^ (sin LF sin £) 


(223) 


where 

5 = 5 - o o 

R' is defined by Equation (220) 


(224) 
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Theoretical burnout velocity and losses . — For trajectory 
and performance optimization studies, it is convenient to know 
the theoretical burnout velocity possible and the velocity 
losses due to gravity, aerodynamic drag, and atmospheric back 
pressure upon the engine nozzle. These quantities may be com- 
puted as follows; 

,Tl heoretical Velocity 

t2 

V C T VAC 

theo - / dt 

J m 

t x '225) 


Speed Loss Due tc Gravity 


V 

grav 



(226) 


Speed Loss Due to Aerodynamic 



(227) 


Speed Loss Due to Atmosphere 
Back Pressure Upon the Engine Nozzle 



Maneuvering Losses 

V - / (c- . - l) « 

t l 


(228) 


(229) 


The resultant velocity Vg(t 2 ) is obtained by adding the components 
computed to the initial value Vg(t^) 

Vg(t 2 ) - V^) + V theo + V grav + V D + Vp (230) 


The maneuvering - losses are valid only if At is zero for the 
engine. 
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Orbital variables and satellite target . — Certain functions 
of use in orbital trajectory calculations have been added to the 
point mass equations of motion used in the Steepest Descent Opti- 
mization Program. These functions permit the specification of 
terminal conditions in inertial space when this is convenient. A 
further set of functions will permit rendezvous calculations with 
a satellite in a circular orbit about a central planet. 

Orbital variable calculations commence immediately after 
the calculation of vehicle inertial velocity. Flight path 
angles in inertial space are computed from the expressions 

-i / Y e + “p cos *l \ 

01 ' tan \ x g / < 231 > 

’= • &) 


The inclination angle, i, is the angle between the plane con- 
taining the velocity vector and the center of the earth, and 
the equatorial plane. 


-L 



Figure 17 .- Orbital Plane Geometry 
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Applying spherical trigonometry to Figure 17 , we obtain the re- 
lationship 


cos i ■= cos <f> T sin o (233) 

u I 

The difference in longitude between the vehicle and the ascending 
i. ode, v, is given by 

tan v = S in $ L tan a J (234) 

The inertial longitude is given by 

6 I = 0 L " V (235) 

and the inertial longitude of the ascending node by 

« = 01 - v (236) 

It is convenient to know the central angle, u, in the orbital 

plane . Measuring from the ascending node , we obtain 

tan $ r 
tan u = tL 

cos Oj (237) 


The orbital variable calculation introduces positional and 
velocity information from a second body. This body is a satel- 
lite considered in a circular orbit about the earth. Its orbital 
height, h s , is specified and remains constant. Position in the 
orbit is computed from an initial central angle, (f>s 0 » by the 
expression 

4>s = ^s 0 + “s* (238) 


The satellite angular velocity is 

inertial velocity, V c , where 

s 

v_ = */ 7 — . ~ 

c s * (R e 4h s ) 


obtained from the satellite 


(239) 


where y g is the gravitational potential constant and Re the 
earth radius. Jt should be noted that Equation (239) assumes 
a spherical earth; for the earth radius is taken as constant, 
and none of che higher order gravitational harmonics are in- 
cluded. Knowing V Co , it follows that 

- v c . 

- u s 

R e +h s (240) 

The variables of this section provide sufficient information to 
either rendezvous with or terminate the trajectory in a speci- 
fied position relative to the satellite. 


62 



VEHICLE CHARACTERISTICS 


Methods by which the aerodynamic , propulsive, and physical 
characteristics of a vehicle are introduced into the computer 
program are presented ir ;his section. Form and preparation of 
the input data are discussed, together with methods by which 
stages and staging may be used to increase the effective data 
storage area allotted to a description of the vehicle's proper- 
ties. 


Aerodynamic Coefficients 

The primary objective of the aerodynamic data input sub- 
program is to provide for a complete accounting of the various 
contributions to the aerodynamic forces and moments regardless 
of the flight conditions of the vehicle being considered. Two 
techniques are available for use in the digital computer pro- 
gram: (a) an n-dimensional table look-up and interpolation and 

(b) an m-order polynomial function of n variables prepared by 
"curve fit" techniques. In the first method, the proper value 
for each term is obtained by an interpolation in "n” dimensions 
where the number of dimensions is taken to be the number of 
parameters to be varied independently plus the dependent vari- 
able. This method has the advantage of accurately describing 
most non-linear variations with reasonable preparation effort. 

The amount of storage space which must be allocated to such a 
method, however, can achieve unreasonable proportions and may 
require substantial computing time for the interpolation as the 
number of dimensions are increased. The second method has essen- 
tially the opposite characteristics; that is, a large amount of 
data may be represented with a small amount of storage space, 
and computation time is held to reasonable limits, but the data 
variations which may be represented must be regular. A substan- 
tial amount of effort can be required for the preparation of 
data by a curve-fit technique. Both these methods are very con- 
venient when the amount of data to be handled is moderate , but 
tend to become unmanageable when large amounts of data are 
required. This usually occurs when the program, having several 
degrees of freedom, is committed to one or the other of these 
two techniques. Therefore, the computer program incorporates 
both of the techniques discussed as a compromise to take advan- 
tage of the more desirable features of both. To do this, a 
general set of data equations have been programmed which define 
each of the aerodynamic forces. In general, the coefficients 
for these equations will be obtained from a curve-read inter- 
polation. Several simplifications may be made to the equations 
depending on the flight condition and vehicle to be considered. 
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Often the particular application will not require some of 
the terms listed in order to describe the flight path and vehicle 
under consideration. The subprogram is arranged so that the com- 
puter will assign a constant value to any curve for which the 
data has not been supplied. For most curves, the constant value 
will be zero. This technique may be used to reduce the time 
required for the preparation of data. Values intermediate to 
those introduced in a tabular listing will be obtained by linear 
interpolation . 


Aerodynamic Forces 

Aerodynamic forces are customarily defined by three mutu- 
ally perpendicular forces. These are lift (L) , drag (D) , and 
side force (Y) . Lift force is perpendicular to the velocity 
vector in a vertical plane; drag force is measured along the 
velocity vector but in opposite direction; side force is meas- 
ured in the horizontal plane, positive toward the right, pro- 
vided the bank angle is zero. If the bank angle is not zero, 

L and Y will be rotated by -B A about the velocity vector. 
Coordinates are shown in Figure 18. 


D 



L 


t 
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Figure 18. — Aerodynamic Forces - Wind Axes 
These forces may be expressed in the form 


L = q;V,h; SC L (V,h,a,B) 

(241) 

D = q(V,h) SCp(V } h,a,8) 

(242) 

Y = q(V,h) SCy(V,h,a,8) 

(243) 


where q is the dynamic pressure and S a convenient reference 
area. The aerodynamic coefficients Cjj, Cq, and Cy may be ex- 
pressed in terms of the aerodynamic derivatives. 


64 



C L * C L 0 + C L a 0 + c L a 2 “| a | + c L g | e| 

+ c L 6 26 2 + a|e| ( 244 ) 

C D ~ C D 0 + C D a | a | + c D a 2 0,2 + C Dg | S I 

+ C V e2 + C DaB M| 8 I (245) 

C Y = C Yo + Cy a ja|+ C Ya 2 a 2 + C Yg 3 
+ C Yg 2 3|B| + C Yag |a|g 

(246) 

Alternatively, the aerodynamic derivatives may be expressed 
as tabular functions of Mach number (Mjj) , a, and 3, that is, 
functions of the state variables and the control variables. 


On occasion, it may be convenient to measure the aerodynamic 
forces in the body axis coordinate system introduced in a pre- 
ceding section, pages 49 to 51. In this case, normal force, (nf), 
is measured along the -2 axis, side force (y) along the y axis, 
and axial force (a) along the -x axis, as in Figure 19. 



Figure 19. — Aerodynamic Force in Body Axes 

The specification of forces in the body axis system is 
similar to that in the wind axis system 


qSCjj 

(247) 

qSCA 

(248) 

qSCy 

(249) 
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where the body axis aerodynamic coefficients are 
C N * C N 0 + Cn, ♦ c N a 2 a 1° | 

♦ C N B t B I + C "e £ ® 2 * C »«6 “ 1*1 

, , j, (250) 

C A * C A 0 + c A a I a I + c A a 2 a 

* C* b Ib| ♦ C* b 2 B 2 * C An8 |a6| 

c y = c y Q + c y u | a l + c y a 2 “ 2 (251) 

+ c y e e + c yg 2 B I 8 I + 'Vae M* 

(252) 


Thrust and Fuel Flow Data 

The techniques employed to introduce thrust and fuel-flow 
data into the equations of motion are developed in an approach 
similar to that employed for aerodynamic data. An n-dimensional 
tabular listiug and interpolation technique is used with the 
independent variables being defined by the type of propulsion 
unit being considered. For the present formulation, the propul- 
sion units are grouped into the following options: (1) rocket, 

(2) air breathing engine. 

Propulsion opuion (1) rocket . — The thrust of a rocket 
motor is assumed variable with stage time, altitude, and, if the 
rocket is controllable, it will also vary with throttle setting. 
The altitude effect is determined by the exit area of the 
nozzle, A e , and the ambient pressure, P. If the thrust is 
specified for some constant ambient air pressure, the altitude 
correction can be calculated within the subprogiam. If the 
rocket motor is uncontrolled, the vacuum thrust, in pounds , will 
be introduced by a tabular listing as a function of time, in 
seconds, and corrected as follows: 

T = Max [T yac - PA e , 0] (253) 

The propellant consumption rate is specified by a tabular 
listing in slugs per second as a function of time, in seconds, 
for the single engine options, or computed from the thrust and 
tl.e engine specific impulse, Isp» for the multiple engine options. 

If the rocket is controlled, the propellant mass flow rate m£ 
is introduced by a tabular listing as a function of throttle 
setting. The thrust is then specified by a tabular listing as 
a function of mass flow rate. 
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Propulsion option (2) air breathing engines . — An air- 
breathing engine is strongly affected by the environmental 
conditions under which it is operating. Engines which would 
be grouped in this classification are turbojets, ramjets, 
pulse jets, turboprops, and reciprocating machines. The param- 
eters considered significant in the program are 

(a) Altitude (h-ft) 

(b) Mach number (M^) 

(c) Angle of attack (a-degrees) , and 

(d) Throttle setting (N-units defined by problem) 

Both the thrust and fuel flow are functions of these variables. 
In order to accommodate these variables, a five-dimensional 
tegular listing and interpolation are used to obtain both thrust 
and fuel flow. The thrust has no further correction as the 
-ffects of all parameters are assumed included in the interpo- 
lated value. 


Engine perturbation factors . — The engine options include 
provision for two data scaling factors for use in parametric 
studies. These are in the form 


T " e 13 T VAC + e 14 


(254) 


Components of the thrust vector . — The equations used to 
reduce the thrust vector to its components along the body axes 
are 


and 


T cosX T 


(255) 

-T sinA T 

cos<(> T 

(256) 

-T sinA T 

sin<|> T 

(257) 


<f> i> and At are defined and explained in the control variable 
section. 


Reference weight and propellant consumed . — Rate of change 
of vehicle mass, m, is set equal to the negative of the total 
mass flow rate, -m t . m is integrated to give variation of 
vehicle mass, m. The instantaneous mass is used in the compu- 
tation of the body motion. The reference weight is obtained 
by an auxiliary calculation 

W T = m (32 . 174) (258) 
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The propellant consumed is computed as 

m f = m Q - m (259) 

where iUq is a reference mass input equal to the initial vehicle 
mass 


Stages and Staging 

A problem common in missile performance analyses and 
encountered frequently in airplane performance work is that 
of staging or the release of discrete masses from the contin- 
uing airframe. The effect of dropping a booster rocket or 
fuel tanks is often great enough to require that the complete 
set of aerodynamic data be changed. Configuration changes at 
constant weight, such as extending drag brakes or turning on 
afterburners, may also require revising the aerodynamic or 
physical characteristics of the vehicle. Another use of the 
staging technique is possible with the present computer pro- 
gram which does not involve physical changes to the configu- 
ration; this technique may be used to revise the aerodynamic 
descriptors as a function of aerodynamic attitude or Mach number. 
With this use of the stage concept, accurate descriptions of 
the forces acting upon the vehicle may be maintained over wide 
attitude ranges, if required. 
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VEHICLE ENVIRONMENT 


The models for simulating the environment in which a 
vehicle will operate are presented in this section. This 
environment includes the atmosphere properties, wind velocity, 
and the field associated with the planet over which the ve- 
hicle is moving. The shape of the planet and the conversion 
from geodetic to geocentric latitudes are also considered. 

In the discussions which follow, the descriptions of vehicle 
environment pertain to the planet Earth. The environmental 
simulation may be extended to any planet by replacing appro- 
priate constants in the describing equations. 


Atmosphere 

The concept of a model atmosphere was introduced many 
years ago, and over the years several models have been de- 
veloped. Reference 20 outlines the historical background of 
the gradual evolution of the ARDC model. The original (1956) 
ARDC model (Reference 20) was revised to reflect the density 
variation with altitude that was obtained from an analysis 
of artificial satellite orbit data. This revision is the 
widely used 1959 ARDC Model Atmosphere and is the basic option 
in the present program. 

The advantage of a model atmosphere is that it provides 
a common reference upon which performance calculations can 
be based. The model is not intended to be the "final /ord" 
on the properties of the atmosphere for a particular time 
and location. The atmosphere properties are quite variable 
and are affected by many parameters other than altitude. 

At the present time, the "state-of-the-art" is not advanced 
to the point where these parameters can be accounted for; 
it may be several years before the effects of some parameters 
can be evaluated. 

1959 ARDC Model Titmosphere . — The 1959 ARDC Model Atmos- 
phere is specified in layers assuming either isothermal or 
linear temperature lapse-rate sections. This construction 
makes it very convenient to incorporate other atmospheres, 
either from specifications for design purposes or for other 
planets. The relations which mathematically specify the 
1959 ARDC Model Atmosphere are as follows (Reference 21) 
the 1959 ARDC Model Atmosphere is divided into 11 layers as 
noted in the table below. 
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Layer 

Hh -Lower Altitude 

Upper Altitude 


(Geopotential) 

(Geopotential) 


Meters 

Meters 

1 

0 

11,000 

2 

11,000 

25,000 

3 

25,000 

47,000 

4 

47,000 

53,000 

5 

53,000 

79,000 

6 

79,000 

90,000 

7 

90,000 

105,000 

3 

105,000 

160,000 

9 

160,000 

170,000 

10 

170,000 

200,000 

11 

200,000 

700, OCO 


For layers 1, 3, 5, 7, 8, 9, 10, and 11, a linear molecular- 
scale temperature lapse-rate is assumed and the following 
equations are used: 


.3046h 

1 + .pl3h/63567oc 


2.269681 x 10- 


Meters 

(260) 

°R 

(261) 

°R 

(262) 

Lb. ,/ Ft . 2 

(263) 

Slugs/Ft. 3 

(264) 

Ft • /Sec • 
Ft .2/Sec. 

(265) 


(266) 


For the isothermal layers 2, 4, and 6, the following changes 
are made 


P = ^-^(Hgp-Ht,) 

pa 


(267) 

(268) 
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Values of the temperature, pressure, density, and altitude at 
the base of each altitude layer are listed below along with 
the appropriate values K^, k 2 , and K^. 




lost 



Quantity 

1 2 

3 4 

5 

6 

*1 

-.22556913 x 10 - * 0 

.13846580 X 10** 0 

-.15920187 X 10** 0 

*2 

-5.2561222 0 

11.388265 0 

-7.5961765 

0 

*3 

0 .15768852 X 10*3 

0 .12086887 I 10* 3 0 

.20623442 X 10’ 3 

*b 

518.686 389.986 

389.968 506.788 

500.788 

298.188 

Pb 

2116.2170 472.67599 

5 X. 9 T 5 M 8 8 . 515*578 

1.2180383 

2.1082485 X 10* 2 

*b 

2.37692 X 10' j 7.0611078 X 10-'< 7.7643892 X 10' 5 2.8803201 X 10‘ h 1.3947125 X 10"' 

4.1190042 X 10'" 

■b 

0 11000. 

25000. 47000. 

53000. 

79000. 



Ii lltr 




7 8 

9 

10 

11 

*1 

.24145841 X IQ -4 .88626910 X 10** .75434123 X 10' 5 

.35071476 X 10‘5 

.22212914 X 10“^ 

*2 

8.5411986 1.706239T 

3.4164794 

6.8329589 

9.7613698 

k 3 

0 0 

0 

0 

0 

*b 

298.188 406.188 

2386.188 

2566.188 

2636.188 

Fb 

2.1811754 X 10’ 3 1.5564912 X 10"* 7.5604667 X 10"* 

5.8$m64ii X 10* 6 

2.9769746 X 10 " 6 

p b 

4.2614856 X 10*9 2.232^12* x 10* 10 1.8458849 X 10 

1.3387990 X 10 - 12 

6.1150607 X 10“ l 3 

■b 

90000. 105000. 

l 60000 . 

170000. 

200000. 

Values of the appropriate constants to 
perature equation, Equation (262) , are 

be applied in the tern- 
listed below. 

Hgp(Ka) 

A 

B 

C 

0 

0-90 

1. 

0. 

mm 

- 

90-180 

.75951115 

.17416404 

220,000. 

25,000, 

180-1200 

.93578678 

.2739659S 

180,000. 

140,000, 


U. S. Standard Atmosphere, ? ' 5 2 . --The part of the U.S. 
Standard Atmosphere, 1962, beiovO kilometers geometric 
altitude (295,276 ft. altitude) is defined in the same way 
as the 1959 model — by the hydros'. <■ cic equation and a piecewise 
linear variation of temperature with geopotential altitude. 
Equations (260) to (268) are, therefore, applicable with a 
different set of constants. These constants, based on the 
published tabulation of atmosphere properties (Reference 22) 
at the base altitudes, are presented below. The 1962 model 
uses a different set of relationships above 90 kilometers. 

These have not been included. The tables define 1962 model 
properties between sea level and 295,800 ft. geometric altitude. 
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Values of the temperature, pressure, density, and altitude 
at the base of each altitude layer are listed below along with 
the appropriate values of K 1 , t. 2> and K 3 « 


Layer 


Quantity 

1 

2 

3 

4 

K 1 

-.2255877 x lo -1 * 

0 

.48012406 x 10“ 5 

.12199559 x 10~‘ 

k 2 

-.5255871 x 10 1 

0 

.32844801 x 10 2 

.12202470 x 10 2 

Ko 

0 

.1576958 X 10~ 3 

0 

0 

Tb 

518.67 

389.97 

389.97 

413. 10^ 

Pb 

2116.217 

1*72.6812 

114.3431 

17.22518 

p b 

.2377002 x 10" 2 

.7061512 x 10~ 3 

.1708202 x 10~ 3 

.2429209 x 10 -1 * 

Hb 

0 

10999.1+7 1 * 

Layer 

19999.191 

32354.854 

Quantity 

5 

6 

7 

8 

K i 

0 

-.7383899 x 10 -5 

-.1572230 x 10- U 

0 

k 2 

0 

-.1709562 x 10 +2 

-.8602817 x 10 

0 

K 3 

.1262323 x 10" 3 

0 

0 

.1891214 x 10~ 3 

Tb 

487.17 

487.17 

454.668 

325.170 

?b 

2.302550 

1.226346 

.3766873 

.2106440 

p b 

.2753526 x 10" 5 

.1466:37 x 10 -5 

.4826665 x 10' 6 

.3773977 x 10"? 

Hb 

1*7051.501 

52042.023 

61077-348 

79192.936 


Within the altitude range considered, T and T M (Equation (262)) 
are equal . M 

Atmosphere limitations . — The validity of the 1959 ARDC 
model is limited to altitudes below 700 km; although the pro- 
gram is arranged to extrapolate the relationships to greater 
altitudes, if desired. Extrapolation to greater altitudes is 
accomplished by altering the cutoff altitude. 

At an altitude greater than 2.6 x 10 6 feet, no calculations 
are made, and the program sets kinematic viscosity, speed of 
sound, pressure, temperature, and density to zero. At and below 
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sea level the parameters, pressure, temperature, and density 
are set to the values below. Other terms are computed as normal. 

Pressure = 2116.2170 Lb/Ft 2 (269) 

Temperature = 518.688 °R (270) 

Density - 2.37692 x 10 -3 Slugs/Ft 3 (271) 

£ 

At altitudes between 90 kilometers and 2.6 x 10 feet, the 
speed of sound is set to 846.50255, and kinematic viscosity 
is set to 2.3519252 x 10“^ over density. Other terms ar~ 
computed as normal. 

The 1962 model is limited to altitudes below 295,8000 feet 
(90 kilometers) . It is suggested that zero values be returned 
above that altitude. At and below sea level, the sea level 
values should be employed. When the atmosphere constants are 
determined from the published tabulations at the base altitude, 
the calculated values at intermediate altitudes may not agree 
with the tabulated values to the number of significant figures 
in the tables. This has been allowed for in the 1959 model 
by developing coefficients with the necessary extra precision 
to give agreement between the calculated values and published 
tables at all altitudes. The values calculated by the 1962 
model are good to about four significant figures, which should 
be adequate for most purposes. 

Kinematic viscosity and speed of sound lose their physical 
significance at very high altitudes, and are not normally 
defined by model atmospheres above- 90 kilometers. The constant 
values by the 1959 model option were added to provide data 
required by the aerodynamic heating routine. The aerodynamic 
heating calculation should not be used with the 1962 model 
option above 90 kilometers. The constant values of v and V s 
in the 1959 model will give reasonable values of Mach number 
and Reyxiolds number for use in the aerodynamics calculations 
to altitudes somewhat above 90 kilometers, say 350,000 feet, 
v above which constant aerodynamic coefficients should be used. 


Winds Aloft 

The winds-aloft subprogram provides for three separate 
methods of introducing the wind vector: as a function of 

altitude, a function of range, and a function of time. This 
facilitates the investigation of wind effects for the conven- 
tional performance studies. The wind vector is approximated 
by a series of straight line segments for each of the methods 
mentioned above. 
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Four options are used to define the wind vector in the 
computer program. The thret components of the wind vector 
in a geodetic horizon coordinate system can be specified as 
t* uular listings with linear interpolations (curve , *eads) 
in the following options. 

Wind options (0) . — In this option the wind vector is 
zero throughout the problem. This allows the analyst the 
option of evaluating performance without the effects of wind. 
This option causes the winds-aloft computations to be bypassed. 

Wind option (1) . — In this option the components of the 
wind vector are specified as a function of time. Wind speeds 
are specified in feet per second and time in seconds. 

Wind option (2) . — The three components of the wind vector 
are introduced as a function of altitude in this option. Wind 
speed is specified in feet per second and altitude in feet. 

Wind option (3 ) . — In this option the components of the 
wind vector are introduced as a function of range. Wind speed 
is specified in feet per second and ran^e in nautical miles. 

The range utilized in this computation is the great-circle 
range . 


By staging of the wind option, it is possible to switch 
from one method of reading wind data to another during the 
computer run. Care must be exercised in this operation, however, 
~s the switching will introduce sharp-edged gusts if there arc 
sizeable differences in the wind vector from one option to 
another at the time of switching. This effect should be avoided 
except in cases where gust effects are fating studied. 


Gravity 

This section presents the equations necessary for the 
introduction of the gravity components into the equations of 
motion. These components were determined by taking partial 
derivatives of the gravity potential equation. The potential 
equation adopted has been recommended for use in the Six- 
Degree-of-Freedoir. Flight-Path Study computer program by AFCRC. 
Constants for the potential equation were determined from 
References 23, 24, and 25. 

Spherical harmonics are normally used to define the grav- 
ity potential field of the Earth, References 23 through 26. 
Each harmonic term in the potential is due to a deviation of 
the potential from that of a uniform sphere. In the present 
analysis the second-, third-, and fourth-order terms are con- 
sidered. The first-order term, which would account for the 
error introduced by assuming that the mass center of the Ee rth 
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is at the origin of the geocentric coordinate system is assumed 
to be zero. With this assumption 



(272) 


where P2, P3# and P4 are Legendre functions of geocentric lat- 
itude expressed as 

P 2 = 1-3 sin 2 ^ 

P 3 = 3 sin - 5 sin 3 4 L 

Pl» * 3-30 sin 2 ^ + 35 sin* * L (273) 


The gravitational acceleration along any line is the partial 
derivative of U along that line. At this point, it should be 
noted that the three mutually perpendicular directions in the 
spherical coordinate system are identical (other than sign) to 
those in the local-geocentric-horizon coordinate system which 
is defined previously. Therefore, the acceleration in the <p L 
direction is identical to gx a # and the acceleration in the R 
direction is identical to -gf . Or in the equation form: 
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Collecting terms: 
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where 


Pg = sin <J> L cos <j> L 

2 

Pg = cos <J> L (1-5 sin (J> L ) 

P^ = sin <p^ cos <f>^ (-3 + 7 sin 2 <J> L ) (278) 

Equations (276) and (277) are used in the gravity subroutine 
with the following values recommended for the constants: 

U g - 1.407698 x 10 16 ft. 3 /sec. 2 

R e = 20,925,631. ft. 

J - 1623.41 x 10~ 6 

K = 6.37 x 10 -6 (279) 

It should be noted that these constants and equations per- 
tain to the planet Earth; however, it is possible to use these 
same equations for any other planet. For this reason, the 
values of these constants is an input to the program so that 
the applicable constants may be inserted for the planet under 
consideration. Due to limited knowledge cf the gravitational 
fields of other planets, it is probable that zero values would 
be assigned to some of the harmonic coefficients when the pro- 
gram is used for entry studies on other planets. 

The above equations are applicable to a non-rotating planet 
as the centrifugal relieving effects caused by the planet's 
rotation are included in the equations of motion. In addition, 
the effects of local anomalies must be added if it is desired 
to make a weight-to-mass conversion based on a measured weight. 
The program has the options of retaining the first, third, and 
fourth order terms. 
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WEIGHTING MATRICES 


The perturbation constraint. Equation (5) , limits perturbed 
control histories to the neighborhood of a time varying 
"volume” lying in the neighborhood of the nominal trajectory. 

At any time point in the control perturbation predicted by 
steepest-descent analysis is directly proportional to the prod- 
uct of instantaneous weighting matrix (metric tensor) inverse 
and partial derivative magnitudes. This is readily apparent 
in the unconstrained solution; for in this case, Equation (43) 
reduces to 



Now consider any point t' along the .trajectories , Figures 20 and 
21 in the two variable case; at this point the integrand of 
Equation (5) defines an inclined ellipse in the (s ai , s a2 ) 
plane. In the usual case, when off-diagonal weighting matrix 
elements are zero, the principle axes of this ellipse parallel 
the (s ai , s a2 ) axes, and the length of the major and minor 
axes are inversely proportional to the weighting matrix diag- 
onal elements (directly proportional to the inverse diagonal 
elements) . 

Now the steepest-descent analysis presented earlier in 
this report can be applied to the case of time pulses in the 
controls at any time t'. This type of problem reduces the 
analysis to an ordinary parameter optimization problem, 
References 15 and 16 in the variables s ai and s a2 . In this 
case, the steepest-descent direction is clear, for it lies 
along the line joining the ellipse center and the ellipse tan- 
gency point to the local performance contours. Further details 
may be obtained from Reference 15. It is clear, however, in 
the case considered that by suitably choosing the weighting 
matrix elements at t = t ' , any descending direction can be made 
that of steepest-descent. 

In general, this behavior persists in the case of a time- 
distributed control perturbation. A badly chosen weighting 
matrix can either greatly inhibit convergence rate or, in an 
extreme case, cause convergence failure due to the limited 
accuracy of digital computation. Examples illustrating this 
behavior in the parameter optimization case can be found in 
Reference 17 . 

Some weighting matrix options available in the program 
of References 1 and 2 are described below. 
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Multiple Control Variable Optimization 

The most insidious types of convergence failures are 
those in which the payoff function fails to reach the optimal 
value, while at the same time the terminal constraints are 
achieved. This problem is prevalent among optimization prob- 
lems involving multiple control variables when a unit weighting 
matrix is employed. The reason for this behavior becomes 
apparent from consideration of an optimization problem involving 
two control variables, and a2» where the weighting function 
W(t) is taken as the unit matrix and is consistently more 
powerful than <*2. More powerful implies that a small change 
in ai will produce a greater change in the payoff function 
than an equal change in 02 will produce for the type of per- 
turbation of interest. In this situation, the greater control 
variable perturbation will tend to appear in rather than 03. 
The total perturbation in the first control variable over a 
series of steepest-descent steps will, therefore, always tend 
to be greater than the total perturbation in the second control 
variable , provided remains the more powerful of the two 
control variables, no matter how many steps in the descent have 
been taken. Now, the total required perturbation in control 
variables during convergence from the nominal trajectory to 
the optimum trajectory is purely a function of the particular 
problem under consideration and the nominal path chosen. There 
is no reason for supposing the total perturbation required in 
the powerful control variable to be either greater than or less 
than that of the less powerful one. It follows that when the 
steepest descent process is presented with a situation in which 
the converse is true, i.e., the weaker control variable re- 
quires the greatest total perturbation, there will of neces- 
sity be a high risk of false convergence. 


The argument can be made more specific. First create a 
measure of the total perturbation required during convergence 
from the chosen nominal to the optimal solution for each con- 
trol variable. For example, this can be achieved by separately 
integrating the absolute value of the perturbation required 
along the trajectory, i.e.. 



( 2 "' ' 


where 


* j a ioptiaun } " j a inoninal| 


( 282 ) 


Total perturbation achieved by the steepest descent 
method after C descents can be expressed in the form 
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{APi(C)| - M 6aij(t)|dt 

I ■'to >1 


( 283 ) 


Here 6o^ j (t) is the perturbation of the i th control variable 
in the jth descent at time t. Now suppose that the r fc ^ con- 
trol variable perturbations are consistently greater than 
the sth # by some order of magnitude P, so that 

ia rJ (t) - 0(P) «a,j(t) (284) 

Inverting this 

- 0(-P) 5ar rJ (t) (285) 

On substituting Equation (285) into Equation (283) 


dP 8 (C) 


T C 

J to 1 J-l 


f T 

C 

dt < l 

E| 

J *o 

J-i 1 


or 


S T C 

y I **rj(t) | dt 

to J-l 1 


(286) 


'to J- 
Monotonic Descents 


Consider the case in which the successive control vari- 
able perturbations at any instant are monotonic as the number 
of descents increases. From Equation (286) the total change 
in the s th control variable will always be P orders of magni- 
tude less than that in the r control variable, no matter how 
many descents are made. In this case, dispensing with the 
inequality in Equation (286) 

C 


AP S (C) = 0(-P) 


S T c 

21 

to J-l 


5otj.j(t) dt 


(287) 


The same remarks are true of Equation (283) . On substituting 
Equation (283) with i - r into Equation (287) 

ZP 8 (C) - 0(-P) 2P r (C) (288) 

That is, the total change in the s th control variable after C 
descents depends only on the change in the more powerful con- 
trol variable and the ratio of their powers. In such a case, 
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once the constraints are satisfied, the r^h control variable 
will approach its final history with regularly diminishing 
steps. The final history of the rth control variable may will 
be near optimum. The stn control variable history will, of 
necessity, be perturbed by smaller amounts on each successive 
descent during this_period until it finally approaches its 
limiting value of AP S (<») . It follows from Equation (288) that 


AP 8 (co) - 0(-P) AP r (oo) 

(289) 

In general there is assurance that either 


* 


AP r (co ) * 5F r 

(290) 

or 


AP 8 (go ) - AP* 

(291) 


If the original total perturbation required in the s*h 
control variable, AP| , is P orders of magnitude less than 
that required in the rth control variable, as it might be if 
previous knowledge of the optimum history of the control 
variable were absent. Convergence would tend to be at least 
one order of magnitude worse in the weaker (s^h) control 
variable than the rth control variable. 

If, on the other hand, the total perturbation required 
in the sth control variable had been Q orders of magnitude 
greater than that of the r*-* 1 control variable, we would have 

Ap£ - 0(Q) Al? (292) 

Combining with Equations (289) and (290) ,in this case 

AP s (ao) - 0(-P) Zp£ * 0(-(P+Q)) AP£ (293) 

If the mean perturbation obtained in the st* 1 control 
variable history after the descent is Aot s , and that required 
for convergence is Aa|, then 

TT" Aa g 

Now problems in which the control variable powers are 
in a ratio of 10^:1 are not uncommon in trajectory optimi- 
zation. It is also fairly common to create a nominal trajec- 
tory in which the weaker control variable has, say, ten times 
greater total required perturbation than the more powerful one. 
In such a case, from Equation (294), when convergence is com- 
pleted, the weaker control variable may be practically unper- 
turbed from the nominal history. 
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In practice, the successive descents need not be monotonic. 
It is, therefore possible for the weaker control variable to 
increase its total perturbation while the more powerful control 
variable oscillates. However, it seems reasonable to assume 
that the descent is "almost monotonic." In this sense, the 
above analysis is "almost correct," and, hence, provides at 
least a qualitative insight into the general behavior of the 
steepest descent process with multiple control variables. It 
should also be noted that the arguments of this section hinge 
on the persistence of unequal control variable powers. 

The possibility of failing to converge to the desired end 
constraints is somewhat more remote than that of failing to 
converge the payoff function. The dominant control variables 
for the payoff function are very often the dominant control 
variables for the constraints and, hence, will continue to be 
perturbed until the constraints are achieved. In addition, 
the control variables usually need not be optimized to achieve 
the end constraints. In any case, failure to achieve the end 
constraints is immediately obvious, whereas the only reliable 
method of checking the payoff function convergence is to obtain 
the same result from as different and widely removed as nomi- 
nal as possible or to apply a time dependent equivalent of the 
topologically invariant warping of Reference 15. 


Control Variable Power 

Previously in this section the concept of control variable 
power has been used; specifically this is a measure of the 
ability of a control variable to influence the final value of 
the payoff function. 

It may be recalled from a previous section that the change 
in payoff function is given by 

S T 

MM W 4t ♦ MM {* x(t o ) } (25a) 

to 

When considering changes in the control variables alone, the 
second term can be ignored. 

Suppose at time t' we create a pulse, i.e., a Dirac Delta 
function of unit magnitude in each of the control variables. 

The change in $ produced by these pulses will be 

«(d*) - M<M MW (295) 

where {1} is a unit column matrix. The elements of the row 
matrix X^G indicate the effect of separate pulses in each 
of the control variables. These elements will be referred to 
as the instantaneous payoff function sensitivities, , or 
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control impulse response functions with respect to performance 


KH - iv*>J [ o <*>] 


(296) 


These quantities measure the ability of a control variable to 
affect the terminal payoff function when constraint changes are 
not considered. The instantaneous payoff function sensitivities, 
X^j(t) , are intimately connected with the optimum control 
variable perturbations. In the case of no terminal constraints 
from Equation (43) the optimum control variable perturbation 
is given by 



(297) 


Substituting Equation (296) into Equation (297) 



(298) 


That is, the optimum perturbation varies directly with the 
instantaneous sensitivities and the inverse weighting matrix. 
If the problem being investigated involves terminal con- 
straints, it follows from Equation (43) and (296) that 


M* * MiK} - [°]lv] M * 1 M} 



M -1 M 


(299) 


These results suggest an approach to the problem of false 
convergence. Since the problem is due to small perturbations 
in weak control variables, an inverse weighting matrix based 
on the control variable sensitivities may be employed to accen- 
tuate the weaker control variables. Effectively changing 
the basis of optimization from that perturbation having the 
greatest change in <J> for a given total perturbation magnitude 
to that perturbation having the greatest change in <f> assuming 
all control variables are equally important and must, there- 
fore, be perturbed by a reasonable amount. 


When concerned with terminal constraint variations, the 
above definition of control variable power may be modified. 
In this case interest lies in control variable perturbations 


83 



which improve the payoff function while providing a prescribed 
change in the constraints. These control variable perturbations 
may be either one of two components or "modes." The first mode 
considered is one in which constraints undergo prescribed 
changes with the minimum control variable perturbation possible. 
From Equation (299) this is seen to be when 

(DPjf - (300) 

with a corresponding mode shape of 

Kl ■ v]' x M (301) 


The second mode considered is one in which the payoff function 
is improved while holding the terminal constraints constant. 
From Equation (299) this mode is given by 


- H'M 

/ JK 


V 


- L X **J C 1 **]" I 

which may be written in the form 


(302) 






l 1 **] 1 I 1 **! 


(303) 


where 


K}’[°]’l> a ] 


(304) 


Substituting pulse variations of this second type into Equation 
(25a) and using Equation (296) 


«(d+) « ; [ S *J [»]-ljK} 


DP 2 . 


-L X ^J C 1 **] { x 4+i (305) 




r 



IW 


(306) 


where s<L , the mixed control variable instantaneous sensi- 
tivities, are defined by 




(30V, 
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Substituting Equation (307) into Equation (303) we obtain 









(3 OF) 


From Equation (308) it follows that when the payoff 
function is improved while leaving the terminal constraints 
unaltered, the control variable perturbations at any point 
will vary directly as the product of the inverse weighting 
matrix and the mixed sensitivity matrix. 


Equations (298) and (308) enable establishment of 
rational methods for choosing weighting functions to insure 
payoff function convergence. When limited to diagonal 
weighting matrices to insure reasonable perturbations in 
all the control variable? , one need only increase those 
diagonal elements of W“* corresponding to the weaker control 
variable elements or decrease those corresponding to the 
powerful elements. Further, the elements of { s$ } or { s^ } 
can be used to decide in which class a particular control 
variable belongs at any instant. End point convergence could 
be improved by basing on G'Xwjq , Equation (301). To 
date, this has not been attempted. 


By integrating the absolute value of the instantaneous 
sensitivities over the whole trajectory, a measure of the 
total control variable power is obtained. If the terminal 
constraint variations are ignored 



(309) 


J.f the terminal constraint variations are held to zero 

{^Hl 0 i*^i dt l (3io) 

The elements of these column matrices will be referred 
to as the integrated payoff function sensitivities. 

The integrated payoff function sensitivities based on 
Equation (309) should approach zero at the optimum; those 
based on Equation (310) do not, of necessity, approach zero. 
Either form, in its own way, serves to measure the overall 
ability of a control variable to affect the payoff function 
and is, therefore, a measure of the control variable power 
previously defined. 

If there were perfect numerical accuracy in the steepest 
descent process, the control variable histories would continue 
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to be perturbed until such time as all the c ;ntrol variable 
powers as measured by integrated payoff sensitivities based 
on Equation (310) were zero. In practice, this condition is 
practically impossible to achieve; in fact, it is often dif- 
ficult to reduce these control variable powers by more than 
an order of magnitude when the weighting function is absent. 
This, then, is the basic reason for the weighting function 
matrix, for without one a high risk of failure to converge 
the weaker control variables is present unless foreknowledge 
of the required total perturbations AFj.* is available. 

It will generally be impossible to obtain the desired 
total perturbations AP^* directly, for to do this would 
require a knowledge of the optimum control variable history. 
In lieu of this knowledge, make the assumption that the AP l * 
all have the same order of magnitude. Reasonable conver- 
gence can then be assured by choosing weighting * .atrices 
based on this assumption. Several such weighting matrices 
based on the payoff function sensitivities will be described 
in the remainder of this section; to date, only diagonal 
matrices have been utilized in this manner. 


Weighting Functions Based on Integrated Sensitivity 


Choose a diagonal weighting matrix in the form. 



A ii + B ii 


M 

E s 






(311) 


where M is the number of control variables. 


With equally powerful control variables, the unit matrix 
is obtained with = l, for then 


(312) 

where S^ 1 is a typical sensitivity. 

In the case of unequal sensitivities this form of the 
weighting function will insure that we have perturbations of 
similar orders of magnitude in each of the control variables. 
For example, suppose there are 

r control variables with S^ ■ O(R), 
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s control variables with S * = 0(S), and 
t control variables with = 0(T) 


then 



S*i 
M + 1 



Integrating 


( 313 ) 



(314) 


Partitioning the matrix according to the power of the control 
variables and considering orders of magnitude 


1 1 

8* + Es* 1 


(r+l)0(R) + sO(s) + tO(T) 


s* * Es* 
S* * 

1 

rO(R) + (s+1) 0(S) + tO(T) 
rC(R) + sO(S) + (t+1) 0(T) 

M + 1 

(M + l) 


(315) 


where sf, and s| are typical sensitivities of order R,S, 
and T, respectively. 


Suppose that R>>S and T, then 



OJHJL. 

(m+ IT 


(r_+_l) 

r 

r 


(316) 


In the extreme case when there is one control variable of 
0(R) and several of 0(S) and 0(T), after the descent is 
complete 


^(0°) 

1 

2 ] 


1 

l 



1 ) 

* 


1 

AF t (oo) 


1 J 


( 317 ) 
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if the descent is monotonic. 


In the absence of a weighting 
from Equation (308) 

|aa 2 (t)|-|8^(t)| 


function in the same example, 

(318) 


On integrating 
/“T 

\ {|a« 2 (t)| l at ~ {s*} 

Jt o 


0(R) 

0(s) 

0(T) 


(319) 


so that the weaker control variables would be practically un- 
perturbed in each descent provided the assumption of R>>S, T 
is retained. On summing over the entire descant, it follows 
that the total perturbation in the weak control variables 
will be negligible compared to those in the powerful control 
variables. 


Weighting matrices based on the integrated instantaneous 
payoff function sensitivities S<J| act in a similar manner by 
emphasizing the first term of Equation (303), instead of 

the complete expression. It is difficult to arrive at a 
quantitative result similar to that of Equations (317) and 
(303) for this type of weighting matrix. For the present it 
must suffice to mention that several cases of false convergence 
in the weaker control variables have been eliminated by the 
use of this type of weighting function. 


Weighting Function Based on Instantaneous Sensitivity 

Suppose there is a single control variable a^, and that 
the power of this control variable varies drastically along 
the trajectory. In region A of the trajectory, let the power 
of ai be several orders of magnitude greater than in region 
B. The greater perturbation will tend to appear in region A 
and, should the discrepancy in control power persist through- 
out the steepest descent convergence, the greater total per- 
turbation will always occur in region A. However, the total 
perturbation required in regions A and B are functions of the 
nominal control variable histories created of the region of 
interest and the problem at hand. Therefore , once more a 
false convergence can occur. 

To be more specific, let region A be that region in which 
t Q £ t <. t' and region B be that region in which t* <. t <. T. 
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Let the power of the control variable in region A be 0(P) 
greater than that in region B. In the absence of a weighting 
function the perturbation mode which improves the payoff 
function directly is proportional to the mixed sensitivities 
payoff function sensitivity, . Therefore, 

(320) 

where 6a j (t) is the perturbation at any point in the jt* 1 
cycle of the descent. 


Following the approach used previously, a W matrix can 
be used which will tend to equalize these perturbations. For 
example , 



(321) 


where may is the largest value of s^ along the trajec 

tory (in practice use the maximum value of s£, from the 
preceding descent) . In this case ^ 


Let 


s * 

a ip max 



be O(P) and let s^* • 

aip mm 


be 0 (Q) . 


(322) 


At the point of greatest power 


MtW ~(2)(0(P)) 


(323) 


and at the point of weakest power 



Mt) Bln ~ (l + ofj}) 0(Q) - Of*) + °( p ) 

(324) 

If P »Q 

we, there fore, obtain 



5a max ^ p 
rfajsin 

(325) 

Without 

the we ghting function 



jggSE - o(p-q) 

**min 

(326) 


Therefore, the solution will be limited to extremely 
small control variable perturbations in region A, unless a 
weighting function is used to alleviate the discrepancy in 
control variable power in the two regions. 
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Combined Weighting Functions 


In general , there may be several control variables whose 
individual sensitivities vary drastically both with respect 
to the independent variable and with regard to each other at 
any instant. The variation with t v 'e independent variable may 
be modulated by using the inverse W matrix which will 
equalize the total sensitivity, i.e., the sum of the 
ual control variable sensitivities, at each point along the 
trajectory. A time varying term of the form 


A 


ii 


+ 



(327) 


will achieve this effect. 


The difference in sensitivity between the control vari- 
ables at any instant may be equalized by utilizing a term 
similar to Equation (311) with instantaneous sensitivities 
replacing the integrated sensitivities 


c ii + D ii 



(328) 


Combining Equation (327) and (328) and adjusting the 
matrix so that with equally powerful control variables 
throughout the trajectory, the unit matrix is obtained with 
A ii = B ii " c ii = B ii * 1* the inverse weighting matrix becomes 




2 (M + 1) 


'ii 


+ D 


ii 


M 

£ 




(329) 


Equation (329) is the most general weighting function 
obtained to date; it has been. utilized successfully with both 
the .mixed sensitivities, , and with the sensitivities, 

s. . The first term in Equation (329) tends to equalize 
differences in total sensitivity at different instants. The 
second term tends to distribute the total perturbation at an 
instant more equally between all control variables. 
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STEEPEST-DESCENT STEP-SIZE CRITERIA 


Application of the Steepest-Descent Method is reduced to a 
routine computation once an automatic scheme or control system 
for determining step-size has been devised. The present section 
describes a step-size control system which has consistently pro- 
duced convergence in the calculation of several hundred diverse 
atmospheric trajectory optimization problems in the period since 
1963. The control system is an integral part of the basic com- 
puter program described previously in this report. This program 
is limited to trajectories having all stages, except possibly 
the last, terminated by a fixed value of independent variable, 
time. The major principles of the control system apply to the 
more general problem of optimal staging described in the multi- 
stage analysis section of the present report. 

A control system is basically a set of tests, mainly of a 
logical nature, which determine a suitable step-size and, hence, 
a control variable perturbation for each successive iteration. 

It can be seen from Equations (35b) and (43) that a step-size is 
determined when the control variable perturbation magnitude DP^, 
the amount of each end point error to be eliminated, di|^, and 
the initial state variable value perturbations, Sx(t Q ) are speci- 
fied. No further mention of the initial state variable value 
perturbations will be made, as this type of perturbation has 
not been included in the basic optimization program. 

It might be thought *'ith the solution of Equation 43 available 
that the choice of step-size is a trivial problem. This judg 
ment would be false, for the solution Equation 43 merely pro- 
vides the optimal perturbation for a very small step: that is, 

it is a linearized solution only. In converging to an optimal 
trajectory, the analyst seeking to keep computer time expendi- 
ture within reasonable bounds will require the use of as large 
a step as possible. These large steps inevitably incur signi- 
ficant differences between linearized predictions of system 
behavior and the actual non-linear system behavior in the presence 
of large perturbations. This can be simply illustrated by 
considering the behavior of any one of the optimization functions , 
f, the payoff function or any of the terminal constraints. For 
an infinitesimal perturbation of the control variable, the linear 
prediction of the change in the function based on the adjoint 
method and the non-linear behavior, as determined by an actual 
trajectory integration will be identical and 

df (e) = Af (e) (330) 

where df is the linear prediction and Af is the actual non-linear 
change. Both are functions of step-size, which is assumed small 
and denoted by e. 
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As the step-size is increased, the difference between df and 
Af will tend to grow progressively. For a large enough step-size 
the linear prediction may cease to provide a reasonable guide to 
the actual effect of the perturbation on the optimization function 
and may, in fact, even be of the wrong sense. This behavior is 
illustrated in Figure 22 for a payoff function being maximized. 



Considering non-linear variations of the type denoted by the lower 
of the two solid lines in Figure 22, the best step-size to use is 
that resulting in the maximum payoff function change (Point A) pro- 
vided this is the only function to be considered. 

Since the linear prediction is available from the adjoint 
equations but all non-linear changes require integration of the 
trajectory equations, it is desirable to determine the maximum 
change in payoff function with the least computational effort. At 
Point A there is a considerable discrepancy between linear predic- 
tion of and actual change in the optimization function. Limita- 
tion of step size to a point such as B, where the linear and 
non-linear changes are in close agreement, results in signifi- 
cantly smaller step-size. If the payoff function was the only 
optimization function to be considered, point A would represent 
a reasonable upper limit on step-size, for any greater step-size 
would reduce the payoff function change. If step sizes greater 
than that denoted by C were to be taken, the payoff function would 
actually be degraded rather than improved. 

Figure 22 illustrates one of the major problems confronting 
the analyst seeking to define a reasonable step size. Too small 
a step, as at B, will produce an extremely well behaved conver- 
gence, but an excessive expenditure of computer time results. Too 
large a step, one significantly greater than that of A, results i n 
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Figure 23.- Possible 
Optimization Function 
Behavior 


the linear prediction becoming 
meaningless, and convergence, if it 
does occur, tends to be erratic. 

Again, this tends to cause excessive 
computer time expenditures. The 
problem is further complicated by 
the fact that usually, instead of the 
single optimization function of 
Figure 22, the behavior of several 
functions must be considered. For 
example, with a payoff function and 
four constraints, a typical problem 
may encounter the changes illustrated 
:n Figure 23. In such a case it is 
not clear what step-size is reasonable. 
In view of the wide variety of be- 
havior exhibited by the optimization 
functions, the need for a control 
system which will make consistent 
and logical choices of step-size, 
becomes apparent. 

In this section, it is assumed 
that a single parameter can be chosen 
to define the step-size and that for 
small perturbations the predicted 
changes in optimization functions 
will vary linearly with this parameter. 
For example, choose the parameter as 
follows: let nominal values be available 
for DP2 and {d<|/} and let these values be 
denoted by DP2 and {di^ 0 }. Now take a 
parameter k to determine step-size in 
the manner 


DP 2 * k 2 DP 0 2 (331) 

(dif>) = k { dip o } (332) 

If {6a 0 }is the control variable his- 
tory generated by the nominal choice 
of step-size parameter k=l, then from 
Equations (43) and (44) 

{6a k > = k { 6a Q } (333) 

and 

d$ k = k d<P Q (334) 

It follows that the perturbations are 
linear with the parameter k, as was 
desired. 
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Control System Philosophy 


There are two philosophies which may be followed in most 
complex decision-making situations. A person may attempt to 
reach a conclusion directly by asking what is the correct course 
to follow, or indirectly by asking which courses are not to be 
followed. The steepest-descent step-size control system follows 
the latter course. The direct approach may at first sight appear 
the more attractive method; however, it should be borne in mind 
that it is usually easier tc determine which courses of action 
should not be followed than it is to determine the particular 
course of action which should be followed. 

Major problems involved in the design of a step-size control 
system are failure to converge and false convergence. The first 
type of failure is immediately apparent, but the latter may be 
difficult to detect. For example, suppose that a case involves 
a single constraint which, after the first M iterations, has 
effectively met the desired terminal value. If the constraint 
value is not permitted to drift away from the desired value sub- 
sequent perturbations will be small by virtue of problem non- 
linearity. In a severe case, this will result in behavior easily 
mistaken for convergence. On the other hand, by permitting the 
constraints to drift off the desired value by means of an indirect 
test, this difficulty may be avoided; this type of behavior is 
illustrated by Figure 24. 

In view of the above and similar types of phenomena, the 
step-size control system has been constructed as a group of very 
loose tests, in the sense that a set of almost obvious decisions 
as to step-size magnitude lead indirectly to a choice of ~ ^p. 


Basic Control System Principles 

Second Order System Behavior . — Each iteration commences with 
a trial trajectory; this trajectory is distinguished from the 
final trajectory of an iteration in that the computation of the 
partial derivative matrices, F and G, is omitted. The computation 
of these matrices usually requires somewhat more computer time 
than the trajectory integration itself. 

This trial trajectory is defined by a step-size parameter, 
k, where k=0 denotes the previous final trajectory and k=l denotes 
the step-size magnitude used to obtain the trial trajectory. 

On completing the trial trajectory, the non-linearities of 
the payoff function and constraints are computed. These are non- 
dimensional measures of the difference between the actual and 
linear predictions of the change in these functions. The payoff 
function non-linearity is defined below. 
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Constraint value 
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NL 


- d^> 


(335) 


and the constraint non-linearities are defined as 


II 


(336) 


Here d<J> and dip denote the linear predicted change in <f> and ip , and 
A<p and Aip denote the actual change between the previous final 
trajectory and the present trial trajectory. 

From the previous discussion for a reasonable step in any 
of these variables, the corresponding non-linearity must be 
neither too small nor too great. 

On completing the trial trajectory, an approximation to the 
actual non-linear variation of the optimization functions with 
step-size parameter k, can be obtained from Taylcr's Theorem by 
making the assumption that the behavior of each function is para- 
bolic. The three conditions defining each of the parabolic 
variations are as follows : 


M 

o 


(337) 

II 

i 

(338) 

H^)-{ 

d(d^) ) 
dk ( 

(339) 


The last of these equations follows from Equation (330) . Equations 
(337) and (338) define two points on a parabola; Equation (339) 
equates the predicted linear slope at the first point to the para- 
bolic slope at that point. Applying these conditions, the approx- 
imate non-linear variations are obtained as functions of k. 


d*lk 2 + d<*> Q .k 


and 


A*(k) - - a«^ 

{a*(k)} * | 


k 2 + 




(340) 


(341) 


Now, the value of k which will provide a specified non-linearity 
in the payoff or constraint functions can be found. Substituting 
Equations (340) and (341) into Equations (335) and (336) 


M 



(342) 
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and 



( 343 ) 


That is, the desired value of k for each quantity is the desired 
value of its non-linearity divided by its non-linearity on the 
trial trajectory. A reasonable value for the non-linearity de- 
sired can be obtained from the geometry of a parabolic variation. 
Consider any of the parabolic approximations to the optimization 
function f, as shown in Figure 25. 



Figure .5. — Parabolic Variations 


For a curve such as OAB, the maximum gain in the function occurs 
at A, and if OAB is parabolic, the non-linearity is 

f NL DC ° (344) 

Accordingly, a reasonable non-linearity for the payoff function 
allowing for the approximation involved is about. 45. For the 
constraints a more conservative value of .3 can be employed. If 
the curve is of the nature of OE, these values still provide a 
reasonable step-size guide. With these assumptions, the step- 
size parameter which gives the desired non-linearity for each 
function by use of Equations (342) and (343) can be computed. A 
basic principle of the step-size control system is to base the 
step-size on the optimization function having the largest k. If 
all the desired non-linearities were equal, this would be equi- 
valent to controlling step-size with the function exhibiting the 
most linear behavior. Additional trial trajectories are made 
when the resulting k<.5 or >2, due to the increased possibility 
of the parabolic assumption being in error if it is either extra- 
polated or interpolated too far. For example, consider Figure 26 
where an interpolation from a trial value causes a reduction in 
f(k) rather than an increase. Similarly, in Figure 27, an 
extrapolation has the same effect. It may be noted from Equations 
(332) and (334) that the trial trajectory corresponds to a k«l. 
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If the largest computed k is less than .5 then take another trial 
trajectory with k=.5 and repeat the above logic. Similarly, if 
k is greater chan 2, a trial with k=2 is .ndicated; however, 
before making such a trial, the control system proceeds to various 
other tests which may reduce the value of k; these tests will be 
described later. 



Figure 26. — Danger of Parabolic Interpolation 



To summarize these tests: their purpose is merely to assure 

that at least one of the optimization functions is reasonably 
linear, a modest requirement for a reasonable perturbation. The 
use of non-linearity in the above manner is the first basic prin- 
ciple of the control system. 
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The second basic principle is that of correcting constraint 
errors gradually. There are several reasons for eliminating 
^onst'*aint errors by a small amount on each iteration, xather 
than by attempting to eliminate the entire error in the first 
iteration. 


First, > irking with non 1 inear equations, the large 

steps which a..*- .squired to eliminate the entire constraint error 
will frequently lie outside the linear range. Hence, after a set 
of time consuming trials cf decreasing step-size, the analyst will 
ultimately be reduced to the gradual elimination of the errors. 


Second, it should be noticed from Equation (43) that of the 
control variable perturbation magnitude DP^ an amount equal to 
| dif< 1 j Itfuj, | { d^ } is required to provide the desired constraint 

changes. If this portion of DP^ is too large, the payoff function 
Equation (44) will be primarily the rosult of constraint changes 
rather than an inherent improvement in the trajectory character- 
istics. In this case, there is a danger that the optimization 
will degenerate into a mere terminal constraint search. 


Third, it must be noted that it is possible to introduce 
local extremals into a problem by the means of terminal constraints. 
This becomes cl-_ar from an elementary example in the ordinary 
calculus. Consider the problem of maximizing a function z(x,y) 
which has a single optimal value as in Figure 28. 

Now seek extremal values of z(x,y) subject to a constraint 

•jU,y,z) = 0 (345) 

It is clear from diagram tnat in the particular case 
considered, there are two solutions: one at A and one at B, 

the global optimum being that at B. Now consider the solution of 
this problem by the method of steepest descent commencing from 
Point C. If achieving the constraint is the dominating influence 
in choosing a step, the solution will tend to traverse a path of 
the nature CDA, and, hence, the lower extremal solution wi. be 
located. On the other hand, if in initially choosing the step- 
size, one pays little attention to the constraint, the likelihood 
of traversing a path such as CEB and locating the global extremal 
is increased. 

From this discussion, it is apparent u . t there are sound 
reasons for nrt attempting to eliminate the complete end point 
error at each step; accordingly the control system initially 
attempts to remove constraint errors of magnitude 

{dv> * AC^ • {*} (346) 


where AC , 


is a small non-dimensional quantity. 
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After N iterations, if certain requirements are met, the control 
system will be attempting to eliminate an error of 


{dip} = N • uC^ * {ip} = C^{ip} (347) 

provided N • ZC^<1. When N • AC^l, the amount of constraint 
error to be removed is given by 

{dip} = {ip} (348) 



Figure 28* — Local Extremal Introduced bi dons train t Function 
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This is the second basic principle of the step-size control 
system, the gradual removal of constraint errors in order to 
emphasize the payoff function role, in the initial iterations of 
a descent. 

The two principles of this section are not adequate to insure 
convergence. It has been necessary to add many other logical 
tests to the control system; some of these are described below. 


Secondary Tests 

The two principles outlined above are far from sufficient to 
insure convergence to the correct solution. They must be sup- 
plemented by many secondary decisions , mainly of an indirect 
nature. The more important ones will now be listed, not neces- 
sarily in the order in which they occur in the actual control 
system. For details of the control system as programmed, refer- 
ence must be made to the appropriate flow charts in Reference 2. 

Determination of Step-Size Magnitude for First Trial of 
Each Iteration . — The step-size magnitude DPr>^ used in the first 
trial of each iteration except the first is automatically based 
on the values used in preceding iterations. For the first iter- 
ation, an arbitrary value must be specified by the analyst; this 
value should be chosen on the large size; the control system will 
very quickly determine the correct value by making trial trajec- 
tories . 

After each iteration in a. attempt to inhibit any tendency 
to a gradual decrease in Dp2, the control system determines the 
trial value from the value finally used on the previous iteration, 
DP2_^ using the expression 

DP 0 2 = 2DP 2 _ 1 (349) 

provided certain other conditions have been met. These other 
conditions are as follows: 

(a) That the value given by Equation (349) is at least 
great enough to provide the constraint change being 
sought. If it is not, then, 

DP 2 * Minf J^J^I^j 1 1 (350) 

That is, DP is taken as the minimum value which 
will provide the constraint change desired, provided 
that this is not greater than 10“4 times the expected 
trajectory cut-off, unless this in turn is a smaller 
quantity than that given by Equation (349) . The 
reasoning behind this test is as follows : 
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(i) Assume first that the rules for con- 
straint changes, to be described in detail 
later, are functioning correctly, and, 
hence, DP^ should be no less than the 
quantity [_ d ^J L^J ~ 1 { d V > • 

(ii) In some cases, the constraint changes 

may be excessive; in this case as 10“4 t 
should insure a reasonable size of per- 
turbation, the control uses this value 
as an upper limit. 

(iii) If, however, the previous iteration used 

a perturbation magnitude DPA_i greater than 
10 - 4 t, then the control system uses Dp2_^ 
as the upper limit on step-size magni- 
tude. 


(b) The quantity. 


grad? 


I44 - Ll**j[l**] { 1 4 *t>\ 


(351) 


is the gradient of <j> with respect to DP if the 
constraint changes are zero; that is, it is the 
measure of how close any trajectory is to the 
optimal trajectory having the same end-points. 

Now grad <f> is usually the difference of two very 
large numbers, and these numbers are the result 
of lengthy numerical computations. In this 
situation, small numerical errors can lead to 
the difference between positive and negative re- 
sults for the value of grad<J> when a trajectory 
approaches the optimal trajectory for a particular 
set of end-points . As these may not be the 
desired set of end-points and as grad<J> is essen- 
tially a positive quantity (see Equation 44) , it 
must be recognized that negative values of grad<f> 
merely mean that a trajectory is practically the 
optimal one to the current end-points. All that 
remains in such a situation is to perturb the end 
points towards their desired values. This is 
accomplished by setting 

DP 2 - |d*J L*]' 1 {*4 DP ° 2 s H M’ 1 { d 4 (352) 


or by following the logic of 
inequality is not satisfied. 


(a) above if the 


(c) On occasion, an idiosyncrasy in a particular 
trajectory may cause the step-si 2 e to become 
severely reduced; this will usually be accompanied 
by an excessive number of trials. After six or 
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eight trials, depending on the circumstance, 
the control system will force a final trajectory 
to be computed. In the next iteration, the 
magnitude of the control variable perturbation 
DP^ for the first trial trajectory will then be 
computed by the expression 

DPo 2 - y (Df^) (DF§_ 2 ) (353) 

instead of by Equation (349) . This value is 
used in an attempt to maintain a reasonable 
perturbation magnitude should an excessive 
number of trials occur. 


(d) If sense switch 4 on the IBM 7094 computer 

console is depressed, the step-size magnitude 
for the first trial of an iteration is auto- 
matically given by 

DP Q 2 = 10 _4 T (354) 


This feature can be used towards the termination 
of a calculation as convergence appears complete, 
to artificially attempt a large step. It is a 
safety device designed to provide assurance that 
a false convergence has not occurred; for by 
this means, a complete family of step-sizes will 
be attempted as the step-size is reduced to an 
acceptable level. 


(e) Two limits are placed on a trajectory; a time 
limit at 


and an altitude limit at 


(355) 


h 


min 


(356) 


If time is integrated beyond T max , or the alti- 
tude below hmin> then the trajectory integration 
is terminated. Should either of these situations 
occur on a nominal trajectory, the control system 
assumes that the nominal trajectory is unsatis- 
factory and terminates the convergence immediately. 
If these limits are violated on any trajectory 
subsequent to the nominal trajectory computation, 
it is assumed that this is the result of a per- 
turbation being too great, and another trajectory 
is computed with 

k = .5 (357) 
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at this point, the value of DP which caused the 
violation of a trajectory limit is saved ( DI, g ave ) 
and is then used as the nominal value on the s 
first trial of the next iteration; for the vio- 
lation of the trajectory limit may be of a 
transient nature, and a major object of the control 
system is to keep the step-size as large as 
possible. 


Determination of Step-Size Magnitude After First Trial . — 
After the •First trial, the step-size magnitude is basically 
controlled by the step-size parameter k, according to the ex- 
pression, 

DP 2 = k 2 .DP 2 Q (358) 

There is an exception to this rule when the step-size is "bouncing 
"Bouncing" means that either a value of DP^ equal to or smaller 
than one already demonstrated to be too small or a value of DP2 
equal to or greater than one already demonstrated to be too great 
is again predicted. Figure 29 demonstrates one way this phen- 
omena can arise. Here a value of DP^ has been computed from a 
value of the step-size parameter k^ QW ; a trial is made and the 
extrapolation is made and a value of k corresponding to point C 
is computed. If this value of k is beyond the point 


k = 2k 

*high low 


(359) 



Figure 29. — Step-Size Bounce Induced by Parabolic Approximation 



The control system will compute a new trial trajectory corres 
ponding to a step-size of khigh' i .e.. 


k “ k high 


(360) 


On completing the trial, the controlling function takes on the 
value at B. The parabolic interpolation from this point pre- 
dicts a value at L less than k. and without a "bounce test," 
a trial would be taken with i ° v ‘ 


k - k 


low 


(361) 


and a closed loop established. Accordingly, if a situation of 
this nature arises. Equation (360) is overruled , and the step- 
size magnitude is determined by a midpoint search 


DP 2 




pf L 


(362) 


Limits on Dimensional Travel of Payoff Function . — The step- 
size parameter k is determined by the first principle described 
previously, that is control with the most linear of the optimi- 
zation functions. This decision is over-ruled if such a step 
causes the dimensional travel of any of the optimization func- 
tions to become excessive. Constraints are placed on the 
dimensional travel of the payoff function in the following 
manner: 


(a) If the problem at hand is one involving maximization 
of the payoff function and 


A* >-*, 


ADV 


* 


ADV 


(363) 


(b) If the problem at hand is one involving minimization 
of the payoff function and 


- *ADV “ *ADV 


(364) 


The permissable adverse <p travel magnitude, ((i^, is 
determined by the expression 

■ ■ - (*■& - 

where 4 >n- 1 is the value of the payoff function at the termination 
of the last ite^p.tion, and <t>max i s greatest value of the 
payoff functio absolute value obtained at the termination of 
any of the previous iterations. 

The adverse <f> travel test described above has its basis 
in the principle of emphasizing the payoff function behavior. 
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Problems are often encountered in which, due to the initial ter- 
minal constraint errors, the unconstrained performance, as 
measured by the payoff function, is better on the nominal 
trajectory than it is on the final optimal trajectory. A prob- 
lem of this nature inevitably involves the loss of performance 
during the major portion of the descent. Now the greatest 
obstacles facing the analyst in applying the Steepest-Descent 
method are false convergence and failure to converge in a 
reasonable number of cycles. Both these phenomena are inhibited 
by the adverse <p travel test when performance has to be given 
up in order to achieve the end points; Fig i 30 demonstrates 
how the test inhibits false convergence in a problem of this 
type. Without the adverse <J> travel tests, the first M itera- 
tions are spent in reducing the constraint error at the expense 
of <J>. At that point (A) in the convergence, if all went well, 
emphasis would return to the payoff function and the optimal 
trajectory obtained at point B. This type of behavior is illus- 
trated by the lines OAB. At point B, however, there is a 
risk of false convergence and the descent may continue in the 
manner of OAC. The adverse <J> travel test, on the other hand, 
will not permit the initial rapid loss of <p and convergence with 
this test included is far more likely to be of the nature of the 
broken line OD. 


Again, in a problem where the performance must tend to 
deteriorate as the constraints improve, a very irregular con- 
vergence may result. This is demonstrated in Figure 31; ini- 
tially, a decline in performance occurs as the constraints are 
improved until the point A^ is reached, /it this point emphasis 
returns to the payoff function; a set of seeps which improve 
performance at the expense of the constraint are undertaken 
until the point B^ is reached. Here emphasis returns to the 
constraint and the process repeats. The resulting convergence 
tends to have the appearance of the lines OA^BnAol^ • . ; the 
adverse <J> travel .test inhibits this irregular behavior and 
tends to lead to a convergence of the nature of OC. 

It should be noted that without the second part of the 
decision of Equation (365) , <f> would be unable to change sign; 
if the end points were attainable with <p = 0 , a false convergence 
such as OD would result. This provides a simple example of how 
an over- restrictive rule in the control system can lead to false 
convergence . 


Whenever the <J> travel fails to satisfy the appropriate 
inequality of Equations (363) or (364) , the parabolic assumption 
is applied to compute a value of k that will result in a satis- 
factory step by solving the equation 


resulting in 

k 


(A4>- d4>)k 2 + d<£ .k * ^ADV 

-d$ ± f (d4>) 2 + 4(A4> - d<t>).$ADV 


(366) 

(367) 


2(A</> - d<*>) 
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Constraint Function « Payoff Function 
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The solution must have one positive and one negative root, 
provided the calculation is performed only when the adverse 
travel is too great, as can be seen from Figure 32. When k has 
been computed from Equation (367) it is multiplied by a factor 
of .9 in view of the approximations involved so that finally 
the acceptable value of k, based on adverse <P travel, is given 

■ .*>5 , (368) 


) 


Limits on Dimensional Travel of Constraints . — Rules which 
specify the amount o£ end-point error to be eliminated on each 
trajectory have been given previously. Due to the non-linear 
nature of the trajectory equations and the necessity of attemp- 
ting to take large steps at each iteration, the actual constraint 
changes may differ considerably from those askec for. Accordingly 
another set of rules which specify acceptable limits on the end- 
point travel must be used; it has proven convenient to state 
these rules in the form 

*BW>i d *i $ Z d *L > (3 6 9) 

(370) 


^ d *i 2 - *fwd. >+i-° 


The permissable non-dimensional limits on adverse con- 
straint travel, and favorable constraint travel, 

are functions of the amount of non-dimensional constraint 
error being eliminated, the number of iterations completed, and 
the number of iterations since the particular constraint error 
changed sign. If less than ten iterations have been completed 
since the constraint error changed sign then: 


^FVD * 3 , 


C^, < 


(371) 


* *5 

* 701 ) * 2 * 5 > 


.5 < (ty< l 


(372) 


*BVD “ * 02 5 

'^FWD * 1,5 > 


J 


(374) 


If more than ten iterations have elapsed since the con- 
straint error changed sign, it is assumed that some difrio_xty 
in meeting the constraint exists. In this case, ’J'fwd and ^BWD 
are based on the number of completed iterations. 
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Payoff function change 



Figure 32. — Application of Parabolic Approximation 
to Adverse <p Travel 
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*FVD. 

cWD. 

1 

y> 

V FWD. 

1 

iJj 

*BWD. 

1 


2.5 
.5, 

1.5 
.025, 


N < 20 


N > 20 
36 


(374) 


(375) 


It has previously been indicated that normally 

C* -MACf (376) 

There is an exception to this rule. The exception occurs when 
the step-size magnitude on the first trial of an iteration is 
less than the amount required to provide the desired constraint 
change and when grad<)> is positive. When this condition occurs, 
Cy is successively reduced by ACy until the constraint change 
is less than the amount the DP^ is capable of providing. 

With and specified, the control system merely 

checks in which direction each constraint is travelling and 
computes by the now familiar parabolic approximation what values 
of k, if any, will cause each constraint to reach the boundary 
towards which it is travelling. The method is demonstrated in 
Figure 33 for a constraint which must be increased and has moved 
in the correct direction. 



Figure 33. — Application of Parabolic Approximation, 
Constraint Moving in Desired Direction 
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If the trial point is at A, the parabolic approximation must 
behave in the manner of A 1 A 2 . The solution sought is at A, and 
the negative root may be ignored. If the trial is at B, the 
approximate solution behaves in the manner of Bj^ * and the 
point B is sought. If the trial is at C , the curve behaves in 
the manner of C^C 2 / and there is no real value of k which will 
produce a point on the forward boundary; in this case, the limit 
on k is ignored by setting k^p^ = ». 

In Figure 34 a constraint which must be increased is con- 
sidered; here, however, the motion is adverse. A trial such as 
that at point D indicates that the point sought is at D; the 
negative solution may be ignored. Similar sketches to those of 
Figure 33 and 34 may be drawn for a constraint which must be 
decreased. 



Constraint Moving in Wrong Direction 


Let the value of k which places a constraint at the appro- 
priate boundary on its travel be denoted by k^rpy^. If the solu- 
tion is complex, adopt the convention that h^Ty^ = «, so that 

-d V(d ^) 2 + M A*- d^)*« 

2(A* - d*) 


when 


and 


fyrVL *® * + < ty')'/'TVL < 0 

d^ 2 + M A^- d^) > 0 


TVL 


(377) 

(378) 


where tf) mTrr = , if the constraint has moved in the wrong 

direction 

= 'J't-Tjr.dij) , if the constraint has moved in the right 
direction 
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In Equation (377) the smallest positive root must be taken. 

Conditions for Ignoring Dimensional Constraint Change Test . 

In some circumstances, the limits on end-point travel are ignored. 
For example, if a constraint has been obtained within the accep- 
table limits, ’I'toL' which are specified by the analyst for the 
particular problem, then its end-point travel ceases to be 
monitored unless the constraint » once again drifts outside the 
acceptable limits. This decision is made in order to avoid 
the possibility of limiting step-size magnitude on the basis of 
a constraint which is essentially met, while considerable errors 
still remain in other constraints or a significant amount of 
performance gain remains. 

The limits on constraint travel are also ignored for a 
constraint error which is being reduced more rapidly than another 
constraint error. This is achieved by creating a measure of 
the end point errors at the termination of the nominal trajec- 
tory. These errors are denoted by i(> ERR . W* 1 ®* 1 , after a number 
of iterations all the constraint errors have been halved , ^err 
is also halved. This process is repeated until the computed 
<|»Err are less than ^tol ; at th i s point 'I'err is set equal to ^tOL* 
Any time a particular end-point error is less than its 

dimensional end point travel will not be tested during the 
following iteration. 

It should be noted that if the <j>TOL are zero, a danger of 
false convergence exists; for if the constraints are essentially 
satisfied before the greatest performance is obtained, a situ- 
ation of the nature of that depicted in Figure 24 exists, and 
the limits on constraint travel may inhibit the development of 
performance . 

It should be noted that whenever the controlling function 
(the one with the greatest k base \ on linearity) is a constraint, 
its end-point travel is always checked, for there is no point 
in controlling with a constraint beyond the permissable limits 
on its travel. If the limits on the travel of a controlling 
constraint cause the step-size to >e less than it would be if 
based on linearity and that constraint is within ^ERR» then an 
attempt to control with the next most linear function is made. 

As the limits on the first controlling function travel can then 
be ignored, it is possible that a larger step will result from 
the use of the second controlling function. The larger of the 
two step-sizes obtained in this manner is then used; if necessary 
this process is repeated with the next most linear function, etc. 

Majority Vote Test. — Only those trajectories on which at 
least half the optimization functions of interest improve will 
be considered satisfactory. The optimization functions of in- 
terest are defined as the constraint functions having errors 
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greater than their respective i^err* and the payoff function 
provided the number of optimization functions of interest is 
odd or zero. 

A trial trajectory which fails to satisfy the majority 
vote test is not permitted to lead to the final trajectory of 
an iteration (valid step) . A valid step which fails to satis- 
fy the test is over-ruled by another valid step. In either 
case, the new trajectory is computed with a step-size based on 
k = .5 


Summary . — The control system essentials have been presented; 
it may be noted that wherever possible, the payoff function 
change is emphasized at the expense of the constraint changes. 

Choice of step-size after the first trial is based on both 
linearity and dimensional changes. A careful examination of 
the various tests will reveal that the step-size parameter is 
basically given by the expression 

- ■ ("-fej v w ^ >tetR ) < 3 ”» 

This value of k must then be checked against the bounce test 
and the majority vote test. If 

• 5 < k < 2.0 (380) 

a final trajectory is computed; otherwise, a further trial 
trajectory at the appropriate limit is computed. 

After a final trajectcry, the majority vote test and the 
adverse travel test must both be satisfied; if they are not, 
then the final trajectory is recomputed with a step-size deter- 
mined by k = .5 or on a computed k,j, TVL 

The control system has shown extreme reliability to date 
provided only that the integration technique employed is ade- 
quate for the problem under consideration. Should a convergence 
failure be encountered by an analyst using the program of 
References 1 and 2, it is strongly recommended that the first 
course of action to follow is a critical examination of the 
integration technique and integration step being employed. 
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